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CN ; ABSTRACT 

o 

Context. While it is currently unclear from a theoretical standpoint which forces and processes dominate the formation of high-mass 

Q stars, and hence determine the mode in which they form, much of the recent observational evidence suggests that massive stars are 
born in a similar manner to their low-mass counterparts. 
, Aims. This paper aims to investigate the hypothesis that the embedded luminous star AFGL2591-VLA 3 (2.3xl0 5 L© at 3.33 kpc) is 

forming according to a scaled-up version of a low-mass star formation scenario. 

Methods. We present multi-configuration Very Large Array 3.6 cm and 7 mm, as well as Combined Array for Research in Millimeter 
Astronomy C 18 and 3 mm continuum observations to investigate the morphology and kinematics of the ionized gas, dust, and 
molecular gas around AFGL2591. We also compare our results to ancillary Gemini North near-IR images, and model the near-IR 
to sub-mm Spectral Energy distribution (SED) and Two Micron All Sky Survey (2MASS) image profiles of AFGL2591 using a 
' Monte-Carlo dust continuum radiative transfer code. 

Results. The observed 3.6 cm images uncover for the first time that the central powering source AFGL 2591-VLA 3 has a compact core 
plus collimated jet morphology, extending 4000 AU eastward from the central source with an opening angle of < 10° at this radius. 
However, at 7 mm VLA3 does not show a jet morphology, but instead compact (< 500 AU) emission, some of which (<0.57mJy of 
■ 2.9 mJy) is estimated to be from dust emission. The spectral index of AFGL 2591-VLA 3 between 3.6 cm and 7 mm was found to be 

between 0.4 and 0.5, similar to that of an ionized wind. If the 3.6 cm emission is modelled as an ionized jet, the jet has almost enough 
momentum to drive the larger-scale flow. However, assuming a shock efficiency of 10%, the momentum rate of the jet is not sufficient 
to ionize itself via only shocks, and thus a significant portion of the emission is instead likely created in a photoionized wind. The 
C 18 emission uncovers dense entrained material in the outflow(s) from these young stars. The main features of the SED and 2MASS 
images of AFGL 2591-VLA 3 are also reproduced by our model dust geometry of a rotationally flattened envelope with and without 
CJ\ \ a disk. 

Conclusions. The above results are consistent with a picture of massive star formation similar to that seen for low-mass protostars. 
However, within its envelope, AFGL 2591-VLA 3 contains at least four other young stars, constituting a small cluster. Therefore it 
appears that AFGL 2591-VLA 3 may be able to source its accreting material from a shared gas reservoir while still exhibiting the 
phenomena expected during the formation of low-mass stars. 

Key words. Radiative transfer - Techniques: interferometric - (Stars:) circumstellar matter - Stars: formation - Stars: massive - 
ISM: jets and outflows 
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1 . Introduction lier times these effects may be reduced. Secondly, while the den- 
sities and temperatures in pristine molecular clouds can easily 

03 ' Does the formation of a massive star differ significantly from that explain the formation of low-mass star-forming cores, without 

of a low-mass star? Arguably, all studies of high mass (>8 M Q ) the input of some stabil ising energy such as an external pres- 

star formation are centred upon this question. There are sev- sure micro-turbulence (pcKee & Tanl l2003h. magn etic fields 

eral possible reasons to expect differences at higher masses, one ( e -g- |Hennebelle eTaU|201Il ICommercon etal.ll2011) or radia- 

being that a massive star is thought to continue to accrete ma- tive heating and outflows (e.g.|&umholzetaL||2012|), these con- 

terial after reaching th e zero-age main sequence (ZAMS, e.g. ditions are not conducive to creating a "monolithic" core that 

IHosokawa et aD[2010h . which is a consequence of its Kelvin- does not fragment, from which the forming massive star can ac- 

Helmholtz contraction time-scale being shorter than its accre- crete all of its mass. Instead, massive stars may start off em- 

tion time-sca le. Therefore , several processes such as rad iation bedded in smaller cores that source most of their mass from the 

pressure (e.g. lYol^l2(502h and ionization (e.g. [^toll2002h may surrounding cluster-forming clump (|Bonnell etalj|2011t |Myers| 

halt, decrease or alter accretion on to t he star. However , in the |2011|), implying that massive stars can only form in clusters, not 

earlier stages of protostellar evolution, IHosokawa et al.l (l2010h m isolation. 

also find that the central accreting stars in their simulations be- Placing these theoretical concerns aside initially, and work- 
come bloated due to accretion, so that the effective temperature ing under the hypothesis that massive stars form as a scaled-up 
and UV luminosity of the protostar remains low. Hence at ear- version of low-mass star formation, in this paper we aim to probe 
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the circumstellar environment of the massive star-forming re- 
gion AFGL 2591 by combining multi- wavelength observations 
and modelling to determine whether any of the signatures of its 
formation differ significantly from those of low-mass protostars. 

AFGL 2591 is a well-studied example of a luminous star- 
forming region (2.1-2.5xl0 5 L^ at 3 33 kpc, lLada et al.l Il984t 
iHenning et aDll990t iRvgl et alJl2012b . One of its most promi- 
nent features is a o ne sided conical reflection nebula observed 
in the near-IR (e.g iMinchin et al.l 11991b iTamura & Yamas hita 
Il992l) . Projected within the Cygnus-X star-forming complex, 
the distance to the source has recently been determined by 
trigonometric parallax measurements to be 3.33+0.1 1 kpc, more 
distant than previously assumed (betwee n 1 a nd 2 kpc, e.g. 
Poetzel et al] 1 19921: iHaseeawa & Mitchell 1 19951: ITrinidad et al.1 



2003; van der Tak & Menten 20Q5)). As will be described below, 



AFGL 2591 actually consists of several objects. However, as one 
source, AFGL 2591-VLA 3, dominates the SED and infrared im- 
ages and hence the luminosity, the name AFGL 2591 will there- 
fore also be used henceforth to refer to this dominant source. 

AFGL 2591 has been studied and modelled by many au- 
thors. For example, one-dimensional modelling of the cir- 
cumstellar g eometry via t h e obs er ved SED has been car - 
ried out by iGuertler etafl (1 19911): Ivan der Tak et al.l (1 19991) : 



Muell er et al.l (l2002bl) an d Ide Wit et al] (12009). Improving on 
these, iPreibisch et all §003 ) used two models - one of a disk 
and the other of an envelope with outflow cavities - to reproduce 
the 40 AU diameter (assuming d= 1 kp c) bright disk of em ission 
observed in their K band image, and lTrinidad et al] (120031) have 
modelled the millimetre emission from AFGL 2591-VLA 3 as 
an optically thick dis k without an env e lope. As part of their 
comprehensive study, Ivan derTaketal.1 d 1999b modelled their 
observed molecular lines (CS, HCN, HCO + ) including an out- 
flow cavity in a power-law envelope, finding that a half-opening 
angle of 30° was able to better reproduce the line profiles, 
van derWiel et al] (1201 lb have modelled a set of six lines de- 
tected toward AFGL 2591 as part of the JCMT Spectral Legacy 
Survey, finding evidence of a cavity or inhomogeneity in the en- 
velope on scales of < 10 4 AU. In addition, as well as observing 
a bi-conical outflow structure in 12 CO (2-1) on s cales of 1-2" or 
several thousand AU. lJimenez-Serra et al] (12012b also uncovered 
evidence for chemical segregation within the inner 3000 AU of 
the AFGL 2591-VLA 3 envelope (assuming d=3kpc, similar to 
our assumed distance of 3.33 kpc). 

The ionized gas emission in the region s urrounding 
AFGL 2 5 91 has prev io usly be en observed by [Campbell 
(Il984bh: llbfani et al] d!995h : iTrinidad et all (120031) and 
Ivan der Tak & Mentenl (120051) from 5 to 43 GHz. These ob- 
servations showed that AFGL 2591 is in fact not an isolated 
forming star, uncovering four continuum sources in the region. 
The observed fluxes of two of these, VLA 1 and VLA 2, gave 
spectral indices consistent with optically thin free-free emission 
from HII regions, and a third source VLA 3, was measured 
to have a steeper spectral index, possibly indicating optically 
thick emission. VLA 3 is also coincident with the central illu- 
minating source of A FGL 2591 observed at shorter wavelengths 
(ITrinidad et al.ll2003|). A four t h radio contin uum source has also 
been detected bv ICampbeiil (Il984ah and iTofani et al] (1 19951) 
(their source 4 and n4 respectively), which shall be referred 
to as VLA 4 in the following sections. In addition, the 3.6 cm 
images presented in this work uncover a fifth source, VLA 5, to 
the south west. 

Knots of H2 and [SII] H erbig Haro objects have 
been detected towa rd AFGL 2591 (ITamura & Yamashitdfl992l: 
iPoetzel et al] 1 1992b . suggesting the presence of shocked gas. 



These are coinc i dent with an east- west bip olar outflow (e.g. 
lLada et al] 1 19841: lHasegawa & Mitchelll 119951) . which extends 
across 5' or 4.8 pc at 3.33 kpc, but also contains a more col- 
limated central small-scale component with an extent of 90" x 
20". 

Evidence for the presence of a disk or rotationally flattened 
material around the central source of AFGL 2591 has been found 
by sev eral authors. The nea r-IR imaging polarimetry observa- 
tions of lMinchin et al] (Il99ll) showed that a disk or toroid of ma- 
terial was needed to appropriately scatter the emission. In addi- 
tion, a large (50" x 80") flattened "disk" of material has been ob- 
served perpendicular to the outflow in observ ations of CS lines 
(I Yamashita et al.l 1 19871). At smaller scales, Ivan der Tak et al] 
( 20061) and lWang et al.l (120121) found evidence for a disk of di- 



ameter 800 AU at 1 kpc (corresponding to 2700 AU at 3.33 kpc), 
which exhibits a systematic velocity gradient in the northeast- 
southwest direc tion. This gradient was also found in the SMA 
observations of [Jimenez- Serra et al] ([2012), which they found 
to be consistent with Keplerian-like rotation around a 40 M© 
star. Finally, bot h OH and water ma s ers have been observed 
towar ds VLA 3 (ITrinidad et"aD 120031: iHutawarakorn & Cohenl 
l2005l: ISanna et al.ll2012l) . Th e Very Large Ar ray (VLAfl 22 GHz 
water maser observations of ITrinidad et al] (2003) uncovered a 
maser cluster, which included a ~0.01" diameter shell-like struc- 
ture on the smallest scales. As well as finding largely consistent 
results, the Very Long Baselin e Array 22 GHz water maser ob- 
servations of lSanna et al] (120121) determined that the maser clus- 
ter is arranged in a v-shape, which coincides with the expected 
location of the outflow walls apparent in the near-IR reflection 
nebula. 

In this paper, we have adopted a multi- wavelength approach 
to further probe the circumstellar environment of AFGL 2591, 
and to address the question of whether it forms in a similar man- 
ner to its low-mass counterparts. Building on previous work, 
the modelling presented in this paper includes the first simul- 
taneous radiative transfer model of the near-IR through sub- 
mm SED, and near-IR images, of AFGL 2591 with a three- 
dimensional axisymmetric geometry. In addition, we present 
new multi-configuration VLA 3.6 cm and 7 mm continuum ob- 
servations that for the first time trace an ionized jet at 3.6 cm, as 
well as 13 CO(1-0), C 18 O(l-0) and 3 mm continuum Combined 
Array for Research in Millimeter- wave Astronomy (CARMA0) 
observations, which trace previously unstudied scales within the 
molecular outflow and envelope, to derive a self-consistent pic- 
ture of AFGL 259 1 . 

Section [2] outlines the new observations carried out in this 
work, describing both the VLA 3.6 cm and 7 mm, as well as 
the CARMA 13 CO, C 18 and 3 mm continuum observations. 
Section [3] presents the archival data used in this paper, includ- 
ing near-IR to sub-mm SED, 2MASS photometry, measured 
2MASS brightness profiles and Gemini North near-IR images of 
AFGL 2591. Section |4l describes the modelling of the near-IR to 
sub-mm SED and near-IR images. Section \5\ presents the results 
for the SED and near-IR image modelling, as well as for the cen- 
timetre and millimetre wavelength datasets. Section presents 



1 Now the Jansky Very Large Array 

2 Support for CARMA construction was derived from the states of 
California, Illinois, and Maryland, the James S. McDonnell Foundation, 
the Gordon and Betty Moore Foundation, the Kenneth T. and Eileen 
L. Norris Foundation, the University of Chicago, the Associates 
of the California Institute of Technology, and the National Science 
Foundation. Ongoing CARMA development and operations are sup- 
ported by the National Science Foundation under a cooperative agree- 
ment, and by the CARMA partner universities. 
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Table 1. Summary of VLA observations 



Wavelength 


Configuration 


Date of 


Program 


No. of 


Time 


Pointing centre 


Phase calibrator 






observation 




antennas 


on-source (hr) 


RA. 


Decl. 


flux density (Jy) 


3.6 cm 


A 


2007 Jul 26 


AJ337 


26 (22) 


1.75 


20 29 24.90 


+40 11 21.00 (J2000) 


1.49 


3.6 cm 


B 


2008 Jan 18 


AJ337 


26 (21) 


1.96 


20 29 24.90 


+40 11 21.00 (J2000) 


1.91 


3.6 cm 


C 


2008 Mar 9 


AJ337 


27 (23) 


1.97 


20 29 24.90 


+40 11 21.00 (J2000) 


2.13 


3.6 cm 


D 


2007 Apr 12 


AJ332 


26 (20) 


0.37 


20 29 24.90 


+40 11 21.00 (J2000) 


1.25 


7 mm 


A 


2002 Mar 25 


AT273 


27 (23) 


1.97 


20 27 35.95 


+40 01 14.90 (B1950) 


3.44 


7 mm 


B 


2008 Jan 18 


AJ337 


26 (24) 


3.43 


20 29 24.90 


+40 11 21.00 (J2000) 


3.15 


7 mm 


C 


2008 Mar 9 


AJ337 


27 (25) 


4.18 


20 29 24.90 


+40 11 21.00 (J2000) 


4.46 


7 mm 


D 


2007 Apr 28 


AJ332 


26 (17) 


0.44 


20 29 24.90 


+40 11 21.00 (J2000) 


1.64 



our discussion, which covers the topics of the jet and outflow of 
VLA 3, and the properties of AFGL 2591 as a cluster. Our con- 
clusions are given in Section [TJ 



well as the synthesised beam and map rms for each part of the 
CARMA observations described in Section [2721 are summarised 
in Table [2 



2. Radio Interferometric Observations 

2.1. VLA 3.6 cm and 7 mm Continuum 

Multi-configuration radio continuum observations at 3.6 cm and 
7 mm were conducted between April 2007 and March 2008 
with the VLA of the National Radio Astronomy Observator)0. 
During this time frame, the VLA continuum mode consisted of 
four 50 MHz bands, two of which were placed at 8.435 GHz, and 
two at 8.485 GHz for 3.6 cm, and similarly two at 43.315 GHz 
and two at 43.365 GHz for 7 mm. At 3.6 cm, observations were 
taken in all four configurations of the VLA (A to D), and at 
7 mm, observations were performed in B through D array and 
combined with A array arch ive data previously published in 
van der Tak & Menten (2005). When combined, these observa- 
tions provided baseline lengths between 35 m and 36.4 km, giv- 
ing information on angular scales from 0.24" to 3' for 3.6cm, 
and 0.05 to 43" for 7 mm. 

For each observation, Table [T] lists the observed wavelength, 
configuration, observation date, program code, number of an- 
tennas in the array (with the number of antennas with useful 
data shown in brackets), time on-source, the pointing centre of 
the target, and the flux density of the gain calibrator determined 
from bootstrapping the flux from the primary flux calibrator. The 
number of antennas with useful data was reduced due to anten- 
nas being out of the array for upgrading and testing of new VLA 
capabilities, general hardware issues, and also some residual bad 
data. For all observations, the gain calibrator was 2015+371 and 
primary flux calibrator was 1331+305 (3C286). 

Data reduction and imaging were carried out using the 
Common Astronomy Software Applications (CASAfl package. 
As 1331+305 was slightly resolved, a model was used for flux 
calibration, which has a flux of 5.23 Jy at 3.6 cm and 1.45 Jy at 
7 mm, allowing data at all uv distances to be used. Antenna gain 
curves and an opacity correction (zenith opacity = 0.06) were ap- 
plied for the 7 mm data. The error in absolute flux calibration is 
approximately 1-2% for the 3.6 cm band, and 3-5% for the 7 mm 
band. The data were imaged using a CLEAN multi-scale decon- 
volution routine, using three scales, and Briggs weighting with a 
robust parameter of 0.5. The rms noise in the final combined 
images were 30 and 56yuJybeam _1 for the 3.6 cm and 7 mm 
images respectively, and the synthesised beams were 0.43" by 
0.40", P.A.=43° and 0.11" by 0.11", P.A.=43°. These values, as 



2.2. CARMA 13 CO, C 18 and 3mm Continuum 

CARMA observations at X ~ 3 mm were taken on 10 July 2007 
in D configuration (corresponding to baseline lengths between 
1 1-150 m). A total of 14 antennas were used during the observa- 
tions, five 10-m and nine 6-m in diameter. 

The CARMA correlator was comprised of two side bands, 
placed either side of the chosen LO frequency, which for these 
observations was 108.607 GHz. Both the upper and lower side 
bands contained three spectral windows: one wide and two nar- 
row, giving a total of six windows. The two wide spectral win- 
dows, centred on 106.7 and 110.5 GHz, had a bandwidth of 
500 MHz and a total of 15 channels, and the four narrow spec- 
tral windows had a width of approximately 8 MHz and 63 chan- 
nels, giving a spectral resolution of 122 kHz. The observed lines, 
13 CO(J=1-0) and C 18 O(J=l-0), lay in the two narrow spectral 
windows in the upper side band at 110.201 and 109.782 GHz. 
At these frequencies, the spectral resolution was -0.33 km s" 1 . 

The phase centre for AFGL 2591 was 20 h 29 m 24 s .7 
+40° ll r l 8". 9 (J2000), and the total time on source was 2.5 hr. 
MWC349 was used as a gain and bandpass calibrator, and 
Neptune as a flux calibrator. Data calibration and imaging were 
carried out in CASA. To do this, the data was exported to CASA 
from its original MIRIAD format after applying line length cor- 
rections and Hanning smoothing. Using a flux for Neptune of 
~4.2Jy at 108 GHz (W. Kwon, private communication), fluxes 
of 1.2 and 1.3 Jy for MWC 349 were derived for the upper and 
lower side bands respectively. As there was insufficient signal- 
to-noise on MWC 349 in the narrow spectral windows, band pass 
calibration was only performed on the two wide spectral win- 
dows. However, as the amplitude bandpass shape varies across 
the 8 MHz band by less than 3%, and similarly for the phase (D. 
Friedel, private communication), this should not significantly af- 
fect our results. To image the data, Briggs weighting and a ro- 
bust parameter of 0.5 were used. The beam size and rms noise 
for each observed line and continuum spectral window are given 
in Tabled 



Table 2. Summary of radio interferometric observations 



3 The National Radio Astronomy Observatory is a facility of the 
National Science Foundation operated under cooperative agreement by 
Associated Universities, Inc. 

4 http://casa.nrao.edu 



A 


Array 


line/ 


Synthesised beam 


Map rms 






cont. 


Size (") 


P A. (°) 


(mJy beam -1 ) 


3.6 cm 


VLA-A to D 


cont. 


0.43x0.40 


43 


0.030 


7 mm 


VLA-A to D 


cont. 


0.11x0.11 


43 


0.056 


3 mm 


CARMA-D 


13 CO 


4.4x3.7 


96 


100 


3 mm 


CARMA-D 


c 18 o 


4.5x3.6 


93 


100 


2.8 mm 


CARMA-D 


cont. 


4.9x4.1 


-4.8 


2.2 


2.7 mm 


CARMA-D 


cont. 


4.3x3.5 


95 


2.1 
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3. Archival Data 

3.1. SED 

The observed near-IR to sub-millimetre SED of AFGL 2591 
is shown in Fig. [T] The SED data points were collected from 
the literature; Table [3] lists the wavelengths, flux densities 
and references for the data displayed in Fig. [TJ The 30 data 
points shown between 3.5 and 45/mi were sampled uniformly 
in log-space from a hi ghly p rocessed ISO-SWS spectra of 
AFGL 2591 dSloan et al.ll2003li taken on 7 Nov 1996, observa- 
tion ID 35700734. The uncertainties for the averaged ISO-SWS 
spectrum were assumed to range from 4 to 22% as per the ISO 
handbook volume V0. As mentioned in Section Q] the SED and 
therefore flux densities of AFGL 2591 given in Table [3] were as- 
sumed to be dominated by one source. Eviden ce to support this 
assumption will be presented in Section [5". 2. 21 via a comparison 
of the SEDs of the resolved sources in the region. 



10-' 




A (/im) 

Fig.l. The SED of AFGL 2591, collated from the literature. 
The best-fitting models for an envelope with and without a disk 
(overplotted solid blue and dashed red lines respectively) are dis- 
cussed in Section [5Jl The errors shown are those reset to 10% 
if the error in the measured flux was <10%. 



3.2. Near-IR 2MASS Images 

A three-colour 2MASS dSkrutskie et al.ll2006h J, H, and K, band 
image of the near-IR emission observed towards AFGL 2591 
is shown in the middle panel of Fig. [2 where red, green 
and blue are K s , H and J bands respectively. This Figure 
shows the conic al near-IR reflection n e bula ( previously ob- 
served by e.g. [k leinman n & Lebofskyl 119751: iMinchin et al.l 
119911: iTamura &Y amashital ll992l) . which is consistent with be- 
ing a blue- shifted outflow lobe from the central source of 
AFGL 2591. 

To obtain integrated J, H and band fluxes for AFGL 2591, 
irregular aperture photometry was performed on 2MASS Atlas 
images. Photometric uncertainties on the fluxes are <1%, how- 
ever it is likely the fluxes are more uncertain due to the choice 
of the irregular photometry apertures, and foreground or back- 
ground stars being included within them. The uncertainty due to 
contaminating stars is likely dominated by the 2MASS source 

5 http://iso.esac.esa.int/users/handbook/ 



20292393+4011105 at 20 h 29 m 23 s .93 +40°11 , 10 ,, .56 (J2000; 
easily visible as a point source in the Gemini North near-IR im- 
ages presented in Section \5?2l * as the apertures were placed to 
avoid all other point sources visible in the 2MASS images. The 
fluxes of this contaminating source are 0.01 1, 0.035 and 0.105 Jy 
at J, H and bands respectively, with band being an upper 
limit. Therefore, the uncertainty in the measured 2MASS fluxes 
of the AFGL 2591 outflow lobe was estimated to be -10%. 

3.3. 2MASS Brightness Profiles 

Flux profiles of AFGL 2591, summed across strips aligned with 
and perpendicular to the source outflow axis (P.A.=256°, taken 
from ISanna et al.l (120121) . with thicknesses of 24.6" and 104" 
respectively), were measured for the three 2MASS bands. The 
background subtracted, normalised profiles are shown in Fig. [3] 
along with the best-fitting models to both the SED and profiles, 
for the models consisting of an envelope with and without a disk 
(blue solid and red dashed lines respectively), which will be dis- 
cussed further in Section 15.11 The background was determined 
by finding the average value within two strips either side of the 
main profile. The errors in the profiles shown in each panel of 
Fig.[3]reflect the uncertainty due to background fluctuations, and 
are calculated as the standard deviation of the profiles measured 
in the two background strips, which were assumed to contain 
minimal source flux. As the strips used to calculate the uncer- 
tainty in the perpendicular profiles contained several bright stars, 
iterative sigma clipping was performed, i.e. pixels lying more 
than 5 standard deviations away from the mean were then ig- 
nored and the mean and standard deviation were recalculated 
and the process was repeated until no pixels were rejected. 

3.4. Near-IR Gemini North NIRI Images 

High resolution Near-IR images of AFGL 2591 were obtained 
from the Gemini website, which were taken under program GN- 
2001A-SV-20 as part of the commissioning of Gemini North's 
NIRI instrument, and were available under public releas^. The 
total integration time was 2 min for J band and 1 min at H and 
K' bands. We aligned the ima ges using Chandra X-ray sources 
in the field (lEvans et al.ll2010l) . giving a positional accuracy of 
0.6", and calibrated the flux level of the images in MJy sr" 1 us- 
ing aperture photometry of bright sources in both the Gemini 
North and 2MASS images. The FWHMs of the PSFs in the 
images range between 0.3-0.4". The peak position of the cen- 
tral source of AFGL 2591 in the J band image was found to be 
20 h 29 m 24 s .86 +40°11 , 19".5 (J2000), which was taken to be the 
position of the powering object. As these images were saturated 
in K' band, 2MASS images were instead used for the modelling 
described in Section [4] However, in Sections 15.21 and 15.31 these 
high-resolution images are compared to the radio continuum and 
molecular line data observed for this work. 



4. SED and Near-IR Image Modelling 

4.1. The Dust Continuum Radiative Transfer Code 

The SED of AFGL 2591 was modelled between near-IR and 
sub-mm wavelengths using the 3D M onte Carlo dust radiative 
transfer code Hyperion described in Robitaille (1201 ll) and the 



6 http://www4.cadc-ccda.hia-iha.nrc- 
cnrc . gc. ca/gsa/s v/dataS VNIRI.html 
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Table 3. Observed near-IR to sub-millimetre fluxes for AFGL 2591, collated from the literature. 



Wavelength 


Flux Density 


Aperture 


Description / Origin 




Qim) 


(Jy) 


radius 








1.235 


0.11 + 10% 


Irregular, < 40" 


2MASS J band (1) ' 


i 




1.662 


0.43 + 10% 


Irregular, < 40" 


2MASSHband(l) 


a 




2.159 


2.30 + 10% 


Irregular, < 40" 


2MASS K s band (1) a 




3.500 


33.68 + 1.347 


ISO apertures 


Sampled ISO-SWS 


spectrum 


(2) 


3.822 


68.98 + 2.759 


ISO apertures 


Sampled ISO-SWS 


spectrum 


(2) 


4.174 


108.4 + 7.590 


ISO apertures 


Sampled ISO-SWS 


spectrum 


(2) 


4.558 


127.2 + 8.906 


ISO apertures 


Sampled ISO-SWS 


spectrum 


(2) 


4.978 


186.8 + 13.08 


ISO apertures 


Sampled ISO-SWS 


spectrum 


(2) 


5.436 


233.8 + 16.37 


ISO apertures 


Sampled ISO-SWS 


spectrum 


(2) 


5.937 


245.0+ 17.15 


ISO apertures 


Sampled ISO-SWS 


spectrum 


(2) 


6.483 


273.7 + 19.16 


ISO apertures 


Sampled ISO-SWS 


spectrum 


(2) 


7.080 


353.2 + 24.73 


ISO apertures 


Sampled ISO-SWS 


spectrum 


(2) 


7.732 


406.2 + 28.43 


ISO apertures 


Sampled ISO-SWS 


spectrum 


(2) 


8.444 


221.6 + 15.51 


ISO apertures 


Sampled ISO-SWS 


spectrum 


(2) 


9.221 


78.69 + 5.508 


ISO apertures 


Sampled ISO-SWS 


spectrum 


(2) 


10.070 


75.49 + 5.285 


ISO apertures 


Sampled ISO-SWS 


spectrum 


(2) 


10.997 


154.3 + 10.80 


ISO apertures 


Sampled ISO-SWS 


spectrum 


(2) 


12.009 


416.8 + 50.01 


ISO apertures 


Sampled ISO-SWS 


spectrum 


(2) 


13.115 


583.4 + 70.01 


ISO apertures 


Sampled ISO-SWS 


spectrum 


(2) 


14.322 


600.0 + 72.00 


ISO apertures 


Sampled ISO-SWS 


spectrum 


(2) 


15.641 


556.8 + 66.81 


ISO apertures 


Sampled ISO-SWS 


spectrum 


(2) 


17.081 


503.2 + 50.32 


ISO apertures 


Sampled ISO-SWS 


spectrum 


(2) 


18.653 


542.3 + 54.23 


ISO apertures 


Sampled ISO-SWS 


spectrum 


(2) 


20.370 


717.7 + 93.31 


ISO apertures 


Sampled ISO-SWS 


spectrum 


(2) 


22.245 


886.1 + 115.2 


ISO apertures 


Sampled ISO-SWS 


spectrum 


(2) 


24.293 


1118 + 145.3 


ISO apertures 


Sampled ISO-SWS 


spectrum 


(2) 


26.530 


1391 + 180.9 


ISO apertures 


Sampled ISO-SWS 


spectrum 


(2) 


28.972 


1728 + 293.8 


ISO apertures 


Sampled ISO-SWS 


spectrum 


(2) 


31.639 


2118 + 465.9 


ISO apertures 


Sampled ISO-SWS 


spectrum 


(2) 


34.552 


2413 + 530.9 


ISO apertures 


Sampled ISO-SWS 


spectrum 


(2) 


on /-too 

37.733 


2810 + 618.3 


ISO apertures 


Sampled ISO-SWS 


spectrum 


(2) 


41.207 


3384 + 744.5 


ISO apertures 


Sampled ISO-SWS 


spectrum 


(2) 


45.000 


3879 + 853.5 


ISO apertures 


Sampled ISO-SWS 


spectrum 


(2) 


60 


4830 + 290 


6' 


IRAS 60//m (3) b 






100 


4390 + 290 


6' 


IRAS 100/rni (3) b 






350 


440 + 20% 


60" 


CSO (4) 






450 


107 + 50% 


60" 


JCMT (5) c 






850 


15 + 20% 


72" 


JCMT (5) c 







Notes. (a) Image shown in Fig. |2] Fluxes found by irregular aperture photometry. {b) Fluxes found using radial aperture photometry on IRIS 
coadded images. (c) Fluxes found using radial aperture photometry on SCUBA Legacy Catalogue images. 

References, (given in parenthese s in the column Descr iptio n / Origin): (1) Skrutskie et al.\ $2006); (2) \Sloan et al\ \200&) ; (3) 
Miville-Deschenes & Lagache (2005); (4) Mueller et al. (2002a); (5) Di Francesco et al. (200$). 



genetic search algori thm used to model IRAS 20126+4104 in 
I Johnston et~aD (l201lh . 

For a given source and surrounding dust geometry sam- 
pled in a spherical-polar grid, Hyperion models the nonisotropic 
scattering and thermal emission by/from dust, calculating r adia- 
tive e quilibrium dust temperatures (using the technique of lLucyl 
1 19991) . and producing spectra and multi- wavelength images. As 
dust and gas are assumed to be coupled in the code we are able 
to probe the bulk of the material which surrounds AFGL 2591. 
In this section, we describe the density structure of the disk, en- 
velope and outflow cavity in the model. 

We model the circumstellar geometry of AFGL 2591 with 
the same 3D axisymmetric envelope with disk geometry success- 
fully employed to model th e SEDs and scattered l ight images 
of low-mass protostars (e .g. iRobitaille et al.l 120071: IWood et all 
2002; Whitnev et al.l2003l) . The axisymmetric three dimensional 
flared disk density is described between inner and outer radii 



R^mdRwk by 

Pdisk(^, z) =Po(^) exp J - ~ 



him) 



(1) 



where 7?o=100AU, tu is the cylindrical radius, z is the height 
above the disk midplane and po is set by the total disk mass 
^disk, by integrating the disk density over z, tu the cylindrical ra- 
dius, and (p the azimuthal angle. The scaleheight h(zu) increases 
with radius: h(rn) = ho (zu/Ro)/ 3 where ho is the scaleheight at 
Ro. We have assumed the parameter = 1.25, derived for irra- 
diated disks around low mass stars (ID' Alessio et al.| [l998). and 
have taken a = 2.25, to provide a surface density of £ ~ tu~ x . 
This equation is esse ntially the same as equation (1) given in 
iJohnston et al.l (120111) . however their second term is removed as 
this pertains to accreting a-disks, whereas our disk is not accret- 
ing. Further, equation (Q} is expressed in terms of Ro instead of 
R+ so that the scaleheight is defined at a radius which is more 
conceptually accessible. 
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Fig. 2. a) Model three-colour J, H and K s band image for the envelope with disk model (RGB: K, H, and J bands respectively), b) 
Observed 2MASS three-colour image of AFGL 2591. c) Model J, H, and K s band image for the envelope without disk model. The 
model images have been normalised to the total integrated fluxes given in Table [3] in order that the morphology of the emission can 
be easily compared. Stretch: red: K s band, 130-300 MJy sr" 1 ; green: H band, 95-150 MJy sr" 1 , blue: J band, 30-60 MJy sr" 1 . 



The density of the circumstellar envelope is take n to be 
that of a rotationally flatte ned collapsing spherical cloud (lUlrichl 
119761: iTerebev etal.ll 1984 with radius tf™ x , 



PenvO,//) =Penv0 



Rr 



-3/2 



1 + 



-1/2 



-l 



(2) 



where p env o is the density scaling of the envelope, r is the spher- 
ical radius, R c is the centrifugal radius, \i is the cosine of the 
polar angle (// = cos 0), and yUo is the cosine of the polar angle of 
a streamline of inf ailing particles in the envelope as r —> oo. The 
equation for the streamline is given by 



l4+fM>(r/R c -l)-fi(r/R c ) = 



(3) 



which can be solved for //o. In our model, we assume the 
centrifugal radius is also the radius at which the disk forms 

Rc - ^disk- 

We chose to describe the envelope density via a density scal- 
ing fact or instead of the envel ope accretion rate and stellar mass 
used in iJohnston et al.l (1201 lb , as this does not require the as- 
sumption of evolutionary tracks - which are currently very un- 
certain for massive stars - to provide a physical stellar radius 
and temperature for a given stellar mass. Instead, we have var- 
ied the stellar radius and temperature as parameters (described 
in Section |4~2l) . 

To reproduce the morphology of the emission in the near-IR 
images, we also include a bipolar cavity in the model geometry 
with density p cav , and a shape described by 



(4) 



where the shape of the cavity is determined by the parameter 
Z? cav which we set to be Z? cav = 1.5, and zo is chosen so that the 
cavity half-opening angle at z = 10, 000 AU is # cav . The cavity 
was included by resetting the envelope density to p cav inside the 
region defined by equation © only where the envelope density 
was initially larger. 

The inner radius of the dust disk and envelope R mm is ex- 
pressed in terms of th e dust destruction r adius R SVL \>, which was 
found empirically by Whit ney et al.l (120041) to be 



^sub - R*(T sub /T±) 



-2.1 



where the temperature at which dust sublimates is assumed to be 
^sub = 1600K, and T+ is the temperature of the star. 

For the g enetic algorithm, whic h is described in detail in 
Section 4.3 of Johnsto netaD (l201ll) . the size of the first gener- 
ation N was set to 1000, and the size of subsequent generations 
M was set to be 200. The code was taken to be converged when 
the^f 2 value of the best fitting model decreased less than 5% in 
20 generations. 



4.2. Input Assumptions: Model Parameters 

A set of plausible ranges for parameters describing the model, 
in which the genetic algorithm searched for the best fit, is given 
in Table |U In addition to this, the stellar radius and temperature 
were also required to lie above the zero-age main sequence de- 
fined as lo g Tzams vs. log (0.9 x /?zams) for the solar metallicity 
models of ISchaller et al.l (Il992l) . where the factor of 0.9 is to en- 
sure that enough models are sampled around the ZAMS, rather 
than strictly above the ZAMS. The stellar temperature T+ and 
a stellar surface gravity of log g = 4 were t hen used to select 
a mo del from the stellar atmosphere grid of ICastelli & Kuruczl 
(2004). To sample the cavity half-opening angle # cav , we first 
determined the relation between the half-opening angle and the 
inclination such that the equation of the cavity, equation Q, pro- 
jected onto the plane of the sky produced the right opening angle 
observed in the Gemini near-IR images, 



= arcsm 



0.826 (sin0 cav )* 



COS Or, 



(6) 



(5) 



When sampling the half-opening angle, the inclination was 
first sampled and equation (6) was then used to determine an 
upper bound for # cav , as radiative transfer effects only work to 
make the cavity look larger than it actually is. For instance, light 
can be scattered into the envelope making it look larger, whereas 
because the density of the cavity is constant, there is no definite 
edge of emission seen within the cavity. The largest half-opening 
angle possible, for an inclination of 90°, was therefore 54° at a 
radius of 10,000 AU. 
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Fig. 3. Black error bars: normalised flux profiles for the three 2MASS bands, aligned with (top) and perpendicular to (bottom) the 
outflow axis. Blue and red lines: the profiles of the best-fitting models to both the observed SED and profiles, for envelope with disk 
and without disk models respectively. Vertical grey dashed lines mark the position of the central source. 

Table 4. Assumed ranges for model parameters as input to the genetic search algorithm. 



Parameter 


Description 


Value/Range 


Sampling 




Stellar radius (R Q ) 


0.6 - 500 


logarithmic 




Stellar temperature (K) 


5000 - 50,000 


logarithmic 




Envelope and disk inner radius (R su b) 


1-100 


logarithmic 


nmax 

env 


Envelope outer radius (AU) 


10 4 - 10 6 


logarithmic 


PenvO 


Envelope density scaling factor (gem -3 ) 


10 -21 - 10 -15 


logarithmic 


^disk or R c 


Disk or centrifugal radius (AU) 


10- 10 5 


logarithmic 


^cav 


Cavity half-opening angle at 10 4 AU (degrees) 


0-90 


linear 


i 


Inclination w. r. t. the line of sight (degrees) 


0-90 


linear 


Pcav 


Cavity ambient density (gem -3 ) 


10 -21 - 10 -16 


logarithmic 


M disk 


Disk mass (M©) 


0.1-50 


logarithmic 


h 


Disk scale-height at radius of 100 AU (AU) 


0.2 - 20 


logarithmic 



4.3. Model Fitting 

As part of the genet ic algorithm used here and described in 
I Johnston et al.l (l201l|), the models were fi t to the data using the 
SED fitting tool of iRobitaille et al.l (120071) . where A v , the exter- 
nal extinction, was a free parameter in the fit. The distance was 
set to be 3.33 kpc, and the visual extinction A v was allowed to 
vary between and 40 magnitudes. 

As described in Robit aiUe et al.1 (120071) . the fitting automati- 
cally takes into account the different circular aperture sizes when 
fitting the models. However, in the case of the ISO fluxes, the 
ISO apertures are not circular, so instead the fluxes were de- 
termined by computing images at the ISO wavelengths, and 
summing up the flux weighted by the slit transmission maps 
(found within the Observer's SWS Interactive Analysis software, 



OSIAd) which give the transmission as a function of location on 
the detector. The ISO slits were rotated to the correct orientation 
for the observation and applied to the model. During the SED fit- 
ting, each data point was weighted by the distance in wavelength 
between its two adjacent data points. This was done to allow the 
algorithm to search for a model which reproduced the shape of 
the entire SED, not only the sections with a high density of data. 

The models were also simultaneously fit to the six flux pro- 
files shown in Fig. O measured from the 2MASS images. The 
model flux profiles were found in the same way as the observed 
profiles, by summing the flux in the model images within 24.6" 
and 104" wide strips aligned with and perpendicular to the out- 
flow axis, and normalising them to the total flux. To produce 



7 http://sws.ster.kuleuven.ac.be/osia/ 
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the convolved model J,H and band images, the model images 
were convolved with the FWHM of the 2MASS point- spread 
functions (PSF) for each band: 3.15", 3.23" and 3.42" for J, 
H and bands respectively, which were determined by fitting 
Gaussian PSFs to stars in the observed images. The model im- 
ages shown in Fig. [2] are high signal-to-noise versions of the im- 
ages used for fitting. 

For both the SED and profiles, the errors were reset to 10% 
before fitting if the uncertainty in the measured flux was <10%. 
This was done to take into account non-measurement errors such 
as variability. The combined reduced^ 2 for each model given as 
input to the genetic algorithm was calculated as 

X combined — % SED % profiles 

(7) 

where X$ed * s me best-fitting reduced x 2 value for the SED 
alone, and;^ 2 rofiles is the overall reduced^ 2 for the profile fits. 

5. Results 

5.1. Results of SED and Image Modelling 

In this paper, we m odel A FGL 2591 as a single source; as we 
will show in Section [7.2.21 the second brightest object in the re- 
gion VLA 1 does not significantly contribute to its infrared flux, 
therefore this is a reasonable approximation. 

The genetic code was run three times - firstly with the model 
parameter ranges given in Table |H as input for all parameters, 
which will be referred to as the 'envelope with disk' model, and 
secondly with the parameter ranges given in Table |4] for all pa- 
rameters except the disk mass and disk scale height at 100 AU, 
Mdi S k and ho, which were set to zero and were therefore not 
treated as model parameters. This second run will be referred to 
as the 'envelope without disk' model. The envelope without disk 
model was run to ascertain whether the SED and images could be 
adequately produced without a disk, i.e. a simpler model which 
had two fewer parameters. For the envelope without disk model, 
^disk will instead be referred to as the centrifugal radius of the 
envelope, R c , as these two parameters are interchangeable in the 
models. The third run had exactly the same setup as the envelope 
with disk model, but was used to determine the repeatability of 
the experiment and any parameter degeneracies, and will be re- 
ferred to as the control model. 

The genetic search algorithm codes were stopped after 28, 29 
and 49 generations for the envelope with disk, without disk and 
control models respectively, when the convergence criterion was 
reached. This corresponds to 6545 models or 2940 CPU hours 
for the envelope with disk model, 6278 models or 3774 CPU 
hours for the envelope without disk model, and 10536 models or 
3960 CPU hours for the control model. 

The resulting SEDs and profiles of the best-fitting models 
for the envelope with and without disk runs are shown against 
the data in Figs. Q] and [3) and the model images are compared to 
the observed images in Fig. [2] The parameters of the best-fitting 
models (i.e. the model in each run with the lowest overall x 2 
value after convergence), are also given in Table [5] 

The minimum reduced x 2 for the best-fitting envelope with 
disk model is 14.9, and the corresponding value for the best- 
fitting without disk model is 18.7, therefore the model with a 
disk provides a marginally better fit to the SED and profiles. 
Comparing the model SEDs to the data shown in Fig.[T] it can be 
seen that they are qualitatively very similar, although the model 
with a disk fits the SED slightly better. This difference is partially 
caused by their disagreement in the mid-IR regime, where the 



without disk model does not follow the data as closely. Figure 
[3] also shows very similar fits to the 2MASS image profiles for 
each model, however the without disk model fits some of the 
profiles marginally better, being more centrally peaked. We note 
that although the reduced x 2 for the best-fitting model with a 
disk is lower, and thus provides evidence that an envelope with 
an embedded disk describes the source better than a model with- 
out a disk, the difference is not large enough to prove this conclu- 
sively. We note that neither model provides a good fit to the J and 
H band profiles along the outflow axis. This is likely to be caused 
by the fact that the emission at these wavelengths is not domi- 
nated by the central source, but instead by scattered light in the 
surrounding medium, which in reality is highly non-uniform, as 
seen in Figs. [2] and H Thus to reproduce the images more accu- 
rately, future modelling will have to take such inhomogeneities 
into account, for exampl e the "loops" in the outflow cavity seen 
in the near-IR images (IPreibisch et al.ll2003b . 

To understand how well the parameters of AFGL 2591 were 
determined, we took the best-fitting models for each run and var- 
ied in turn each parameter across the original parameter ranges, 
while holding the other parameters constant. We then fit the re- 
sultant SED and images for each new model, and calculated their 
overall reduced^ 2 -values. From these we could therefore deter- 
mine the x 2 surfaces in the vicinity of the the best-fitting model 
parameters and hence understand their uncertainties. Figure |4] 
shows such a plot for each varied parameter, for all three runs. 
Firstly, we can see the best-fitting values and shape of the x 2 
surfaces for stellar radius and temperature are different for each 
run, which is due to a degeneracy between these parameters. The 
X 2 surfaces for the stellar luminosity, which is a combination of 
these two parameters, is shown as a panel in Fig. |4j this sur- 
face is very well defined and very similar for all runs, showing 
that fitting the SED allows this parameter to be accurately and 
non-ambiguously determined. In the next panel we show the^ 2 
surface for the parameter orthogonal to the stellar luminosity, 
R\ TV, is flat and thus poorly determined. Therefore the SED or 
image profiles do not provide a handle on either the stellar ra- 
dius and temperature, only the stellar luminosity. Similarly, the 
envelope density scaling factor p env o and or R c are degen- 
erate, as they both determine the envelope density, which can 
be seen from equation (2). The final two panels of Fig.|4l show 
the orthogonal parameter combinations p env ^J 2 and p env Rq 2I?> • 
Here p env R?J 2 is similar between the three runs, and has a sharp 
X 2 surface, while the best-fitting p env 2/3 value changes by an 
order of magnitude. Thus we cannot uniquely determine these 
two parameters from our modelling, but we can determine a pa- 
rameter that combines these which governs the overall envelope 
density and mass, p env R 3 J 2 - 

Looking at the remaining parameters, it can be seen that al- 
though the minima in their x 2 surfaces are not sharp, the shape of 
their x 2 surfaces is similar for the three runs, with R mm tending 
to smaller values, R™* tending to intermediate or large values, 
and the inclination preferring intermediate values. The x 2 sur- 
face for the cavity half-opening angle # cav changes between each 
run, however is generally poorly constrained. It is also impor- 
tant to note that the half-opening angle is somewhat degenerate 
with the inclination, as the inclination of a model determines the 
possibl e upp er limit for the half-opening angle as described in 
Section 14.21 The shape of the x 2 surface for the cavity density 
p C av is variable at higher values, but always has its minimum 
at low values, between 10" 21 and ~ 2 x 10" 20 gem -3 . For the 
two runs which include a disk, the x 2 surfaces of the remain- 
ing parameters Mdisk and ho are almost flat across the sampled 
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Table 5. Parameters of the genetic algorithm best-fitting models 



Parameter Description Value for envelope Value for envelope Value for envelope with 

with disk model without disk model disk control model 





Stellar radius (R ) 


90 


200 


52 




Stellar temperature (K) 


13,000 


8900 


18,000 




Envelope and disk inner radius (R SU b) 


2.9 


3.4 


4.3 


omax 
env 


Envelope outer radius (AU) 


180,000 


420,000 


790,000 


PenvO 


Envelope density scaling factor (g cm -3 ) 


2.4xl0" 19 


l.lxlO" 18 


5.6x10" 


^disk Of 


Disk or centrifugal radius (AU) 


35,000 


9500 


2700 


$cav 


Cavity half-opening angle at 10 4 AU (degrees) 


16 


20 


49 


i 


Inclination w. r. t. the line of sight (degrees) 


59 


46 


54 


Pcav 


Cavity ambient density (g cm -3 ) 


2.1X10" 20 


2.0xl0" 20 


1.0xl0" : 


M disk 


Disk mass (M ) 


14 




16 


ho 


Disk scale-height at radius of 100 AU (AU) 


11 




11 


u 


Stellar luminosity (L©) 


2.3xl0 5 


2.3xl0 5 


2.4xl0 5 


PenvO 


Envelope density parameter (g cm -3 AU 3/2 ) 


1.6xl0" 12 


l.OxlO" 12 


8.0x10" 



ranges, although they slightly prefer values close to the upper 
limit. However, for ho the^ 2 rises steeply toward the upper limit 
of the parameter range after its minimum at 1 1 AU. 

In the remainder of this Section, we compare the results 
of our modelling to the recent results from other studies of 
AFGL 2591. Firstly, our modelling finds values of roughly 2- 
5/? su b for the envelope and disk inner radius R mm . Using equa- 
tion (5), we find = 34-39 AU for the three best-fitting mod- 
els, and hence we expect values of ~ 70 - 200 AU. This 
va lue is in agreement with the radius of the inner cavity adopted 
by I Jimenez Serra et al.l (120121) . w ho in turn took this valu e from 
the disk of emission observed by iPreibisch et al.l (120031) . likely 
to be tracing th e inne r rim of the envelope and disk. At 1 kpc, 
IPreibisch et al.l (|2003|) found the radius of this emission was 
40 AU, which is thus 130 AU at 3.33 kpc. 

We compare the maximum envel ope radius R™? to the re - 
sults of the JCMT line observations of van der Wie l et al ] (1201 11) . 
which were able to trace the largest envelope scales, with the 
maximum radius in their model set to 200,000 AU (scaled to 
3.33 kpc). This is in agreement with our SED modelling results, 
which prefer large envelope radii, on the scale of several hun- 
dreds of thousand AU. 

We find the cavity half-opening angle # cav to be poorly con- 
strained by our modelling, although we determined an upper 
limit for the cavity half-opening angl e of 54° in Section 14.21 
(defined at 10 4 AU). IWang et al] (120121) adopted an opening an- 
gle of 60°, and hence a half-op ening angle of 3 °, in their 
modelling. The models u sed in I Jimenez- Serra et ail (1201 2l) and 
van der Wiel et do not include an outflow cavity. 

We find the inclination to prefer intermediate values, be- 
tween roughly 30-65°. This is in agreem ent with the adopted 
inclination of 30° by IWang et akl (1201 2|), who in turn used 
the value determined by Ivan der Tak et al] (120061) . However, 
Jimenez-Serra et al. (2012) found an inclination of 20° by mod- 
elling the position- velocity diagram of combined emission from 
several lines. Thus, previous results appear to favour an incli- 
nation which is slightly closer to face-on. As the inclinations 
for the envelope with and without disk models were 59 and 46° 
respectively, we therefore adopted an inclination of 50° to deter- 
mine several of the physical properties of the jet and outflow in 
Section [67T1 as a compromise between these two values. 

We found low cavity densities p cav were preferred for 
all runs, with best-fitting values between 10" 21 and ~ 2 x 



10" 2Q g cm" 3 . In IWang et al] (120121) and Ivan der Tak et al.l 
(11999 ). their model cavity densities were set to zero. 

Our modelling found the mass of the disk to be on the order 
of 10-20 M Q , whereas using .5" 203.4 GHz dust continuum ob- 
servations Wan g et all (120121) determined the disk mass to be ~ 1- 
3 M . However, as the disk mass was one of the most poorly de- 
termined parameters by our modelli ng, and as the inter feromet- 
ric dust continuum observations of Wa ng et all (1201 2l) directly 
measure the column density and thus mass of the disk, with- 
out the intervening emiss ion or scattering fr om the envelope, we 
suggest that the result of I Wang et al] (120121) is more robust. 

It was not possible to compare o ur best-fitting valu es for the 
disk scale-height ho to the model of IWang et al] (12012b . as they 
instead used a sin 5 6 function to describe the vertical structure of 
the disk, where is the angle from the polar axis, as opposed to 
our Gaussian profile with z. 

The best-fitting luminosity of AFGL 2591 is 2.3xl0 5 L for 
both the envelope with and without disk models respectively, 
which is in agreement with previous results (2.1-2.5xl0 5 L Q 
at 3.3 3 kpc, iLadaetalJ 11984 iHenning et all 1 19901: iRvgl et al.1 
l2012h . 

The envelope density parameter p env o R 3 J 2 is a proxy for 
the mass of the envelope, which dominates the total mass of 
the enclosed circumstellar material at large radii, as it describes 
the scaling of equation (2). We find the total mass of the enve- 
lope plus disk is 1200M© for the envelope plus disk model, and 
2700 M for the envelope without disk mo del. The total mass of 
gas associated with the region observed bv lMinh & Yamzl (2008) 
is 2xl0 4 M o , which is scaled to 3.33 kpc and within a radius 
of 5.8 pc. Thus although our model traces a sizeable fraction of 
the envelope material within a radius of several p c, it may miss 
more diffuse material on the scales observed by Minh & Yang 
( 2008), and thus could underestimate the true total mass asso- 
ciated with the region. To determine whether the mass distribu- 
tion on smaller scales of our best-fitting models was consistent 
with other results, we comp ared the total mass within a radius 
of 2700 AU to the model of iJimenez-Serra et al.1 (120121) . Within 
this radius, which is half the determined si ze of their envelope 
5400 A U, the total mass of the model of IJimenez-Serra et al] 
(1201 2l) is 1 M . Comparing this to our model we find an en- 
closed mass within 2700 AU of 0.1 M© for the model without 
a disk, and 1 M Q for the model with disk. Therefore in the case 
where there is a disk, the enclosed masses are in agreement. 
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Fig. 4. Plots of the x 2 surface along the axis of each model parameter for each run, determined by varying each parameter while 
holding the other best- fitting parameters constant. Red dashed line: envelope without disk run, light blue line: envelope with disk 
run, and black line: control run. 
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Fig. 5. a) Map of the 3.6 cm continuum emission surrounding 
AFGL 2591. Contours are -3, 3, 4, 5, 7, 10, 15, 20, 30, 40, 
50... 100 x rms noise = 30/Jybeam" 1 . Grey scale: -0.03 to 3.77 
mJy beam -1 (1.2 x peak value). The synthesised beam is shown 
in the bottom left corner: 0.43" x 0.40", P.A. = 43°. b) Close- 
up of the 3.6 cm continuum emission towards VLA 3. Contours, 
grey scale and beam as in panel a). 



5.2. 3.6 cm and 7 mm Continuum 

Figure [5]presents the observed multi-configuration 3.6 cm image 
of AFGL 2591. Panel a) shows the entire region surrou nding the 
source, including the four sources first observed by ICamp bell 
(Il984al) at 6 cm (VLA 1 through 4). In addition, a new, low sur- 
face brightness source was detected: VLA 5, which varies in 
peak flux from 3-5 cr across the extent of the source. With the ad- 
dition of shorter baselines, more extended em ission is recovered 
compared to the 3.6 cm A array images of iTofani et al.l (1 19951) 
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Fig. 6. Map of the 7 mm continuum emission towards 
AFGL 2591. Contours are -4, 4, 7, 10, 15, 20, 25, 30 x rms noise 
= 56/Jy beam -1 . Greyscale: -0.06 to 2.14 mJybeam -1 (1.2 x 
peak value). The inset panel shows a close-up of the 7 mm emis- 
sion from VLA 3. The synthesised beam is shown in the bottom 
left corner of both images: 0.11" x 0.11", P.A. = 43°. 



and iTrinidad et al.l (120031) . Most noticeably, the emission from 
source VLA 3, which is coincident with the position of the cen- 
tral source in the near-IR images, is better-determined; a close- 
up of the emission from this source is shown in panel b) of Fig. [5] 
The 3.6cm emission from VLA 3 is consistent with a compact 
core plus jet morphology. The dominant side of the jet extends to 
the east, with a deconvolved width and length of <0.2" and 1.2" 
(<670AU and 4000 AU at 3.33 kpc), position angle of -100°, 
and an opening angle (at a radius of 4000 AU) of <10° derived 
from its deconvolved width and length. The east jet is resolved, 
ending in a 'knot' at 20 h 29 m 24 s .98 +40°11 , 19 ,, .3 (J2000). There 
is also a smaller corresponding jet to the west, seen as a slight ex- 
tension of the source in this direction. The morphology of VLA 3 
is consistent with a jet- system which is approximately parallel 
to th e larger-scale flow (observed by e.g. lHasegawa & Mitchell 
199 5j), and therefore these observations confirm the proposal of 
Trinidad et al. (2003) that VLA 3 is its powering source. 

Figure [6] shows the 7 mm multi-configuration image of 
the region surrounding AFGL 2591, in which only VLA 1 and 
VLA 3 were detected at this resolution. VLA 3 is compact, with 
a radius of < 500 AU, and is only slightly resolved in this im- 
age, showing no evidence of the eastern jet seen in Fig. [5] at 
3.6 cm, which may partially be due to the 7 mm images having 
al most twice the rms noi se. However, the A array-only image 
of Ivan der Tak & Mentenl (120051) . which has a synthesised beam 
size of 43 by 37 mas, shows that the source is slightly extended 
to the south-west, in a similar direction to the counter-jet seen 
in the 3.6 cm image. This extension is also seen in our multi- 
configuration 7 mm image. Using VLBA water maser observa- 
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Table 6. Measured properties of the observed 3.6 cm (8.4 GHz) and 7 mm (43 GHz) continuum sources 



Source A Peak position Peak flux Integrated flux Deconvolved PA. 

name (cm) (J2000) density (mJy beam -1 ) density (mJy) source size (") (degrees) 

VLA1 JZ 20 29 24.59 +40 11 15.0 3.14 + 0.07 80 + 1.6 4.9 x 4.0 120 

0.7 20 29 24.580 +40 11 14.50 0.46 + 0.06 64 + 3.2 1.8 x 1.6 90 

VLA2 3.6 20 29 24.52 +40 11 20.1 0.45 + 0.03 5.2 + 0.11 2.5 x2.3 110 

0.7 

VLA3 3.6 20 29 24.88 +40 11 19.5 0.72 + 0.03 1.52 + 0.03 1.9 x 0.84 90 

0.7 20 29 24.882 +40 11 19.45 1.8 + 0.11 2.9 + 0.14 0.35 x 0.26 50 

VLA4 3.6 20 29 24.32 +40 11 28.0 0.39 + 0.03 0.99 + 0.02 1.1 x0.78 90 
0.7 

VLA5 3.6 20 29 24.00+40 11 10.3 0.15 + 0.03 10.9 + 0.22 4.4x2.9 
0.7 



tions, IS anna et all (120121) found that this emission most likely 
arises from the outflow cavity walls, which the masers they ob- 
served were also found to trace. 

Table [(^provides the measured positions, peak and integrated 
fluxes, as well as deconvolved sizes for the sources observed in 
the 3.6 cm (8.4 GHz) and 7 mm (43.3 GHz) images. These were 
measured using a custom-made irregular aperture photometry 
program, with which the integrated flux density was measured 
above 1 sigma. Uncertainties in the aperture fluxes were calcu- 
lated to be a combination of the uncertainty due to the image 
noise over the aperture, and the maximum VLA absolute flux 
error, which is 2% and 5% for 3.6 cm and 7 mm respectively. 
Similarly, the peak flux uncertainty was found by combining the 
1 sigma flux density with the VLA absolute flux error. The de- 
convolved source size at the 3 - <x level was measured by taking 
the major axis of the source to be along the direction which it is 
most extended, to within a position angle of 10°. 

Figure[7]compares the 3.6 cm emission from the region to the 
near-IR emission recorded in the Gemini North images. The five 
detected HII regions line-up with several features in the near- 
IR image. Firstly, the peak of VLA 3 is coincident with that of 
the central source of AFGL 2591, at the apex of the one sided 
reflection nebula. VLA 1 is instead anti-correlated with the dif- 
fuse near-IR emission. This dent or cavity in the cloud can be 
seen more cle arly in the K band bi spectrum speckle interferom- 
etry image of iPreibisch et al.l ([2003). VLA 2 als o appears to be 
coinci dent with the close binary discovered by IPreibisch et al.l 
(2003), and VLA 4 with a source in the Gemini North image 
at 20 h 29 m 24 s .3 +40°11'28" (J2000; however this source is too 
faint to be seen in Fig.O. In addition, VLA 5 may be powered by 
the bright source at 20 h 29 m 23 s .96 +40° 1 1'09".25 (J2000; which 
is the same source as 2MASS 20292393+4011105, previously 
mentioned in Section [3T2b . The sources VLA 1 and VLA 3 are 
discussed in detail below. 



5.2.1 . The Central Source of AFGL2591 , VLA 3 

Figure [8] shows the SED in frequency space of the central source 
of AFGL 2591, VLA 3, in black error bars. The fluxes shown 
for VLA 3, which are taken from the literature and this work, 
are listed in Tables [3] (for A <1 mm) and [7] (for A >1 mm). At 
A > 1 mm, fluxes have been listed from the observations with the 
largest available uv coverage. 

A greybody was fitted to the fluxes between 100 and 850 yum, 
giving a temperature of 25 +7 K and a dust emissivity exponent 
P of 2.3+0.6. For these four fitted fluxes, the image photometry 
apertures contained most if not all of the source flux. 




25.20s 24.80s 24.40s 20h29m24.00s 

RAQ2000) 



Fig. 7. Three-colour JHK' Germini-North image of AFGL 2591 
overlaid with the 3.6 cm contours from Fig. Stretch of Gemini 
image: red, K' band: 150 - 2500 MJy sr" 1 ; green, H band: 150 - 
500 MJy sr" 1 ; blue, J band: 40 - 80 MJy sr" 1 . 



It can be seen in Fig. [8] that the fluxes at millimetre wave- 
lengths (87-226 GHz) fall short of the flux expected from the 
fitted greybody. However, as these fluxes are taken from inter- 
ferometric observations, it is likely that they only contribute a 
fractio n of the actual flux. For instance, Ivan der Tak & Mentenl 
( 2005) note that their quoted 226 GHz continuum flux represents 
only 5% of the total flux, due to lack of shorter baselines. An ar- 
row is shown on Fig. [^representing this correction, which agrees 
well with the flux expected at 226 GHz from the fitted greybody. 

Extrapolating this greybody to 7 mm, a flux of 2.3 mJy is 
expected, but the measured flux of VLA 3 at 7 mm is 2.9 mJy. 
However, it is likely that a portion of the 7 mm flux is due to 
ionized gas emission as well as dust emission, and a significant 
fraction of this combined emission is resolved out by the inter- 
ferometer. To estimate the contribution from dust emission at 
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Table 7. Overview of observed fl uxes for VLA 1 to 
VL A5. References a r e: (1) ICampbelJ (Il984ah. ( 2 ) this work 
(3) Ivan der Tak et al.1 (119991) . (4) IjVIarengo et al.1 (l2000b . (5) 
iTrinidad etal.1 (l2Q03h 



Source 


Wavelength 


Integrated flux 


Ref. 


name 




density (mJy) 




VLA1 


6 cm 


79 + 2.0 


(1) 




3.6 cm 


80+1.6 


(2) 




7 mm 


64 + 3.2 


(2) 




3.4 mm 


87 + 1.4 


(3) 




2.83 mm 


71 + 1.2 


(3) 




2.81mm 


76 + 15 (20%) 


(2) 




2.7 mm 


72 + 14 (20%) 


(2) 




2.6 mm 


93 + 2.5 


(3) 




18.0//m 


7.96 + 0.8 x 10 4 


(4) 




12.5 /im 


1.59 + 0.2 x 10 4 


(4) 




11.7 /mi 


1.84 + 0.2 x 10 4 


(4) 


VLA 2 


6 cm 


3.61 + 0.72 


(1) 




3.6 cm 


5.2 + 0.11 


(2) 


VLA 3 


6 cm 


0.4 + 0.1 


(1) 




3.6 cm 


1.52 + 0.03 


(2) 




1.3 cm 


1.57 + 0.4 


(5) 




7 mm 


2.9 + 0.14 


(2) 




3.4 mm 


29.5 + 0.8 


(3) 




2.83 mm 


38.7 + 0.7 


(3) 




2.81mm 


71 + 14 (20%) 


(2) 




2.7 mm 


79 + 16 (20%) 


(2) 




2.6 mm 


52.9 + 1.5 


(3) 




1.3 mm 


-151 +4.5 


(3) 




18.0//m 


7.53 + 0.75 x 10 5 


(4) 




12.5/mi 


7.45 + 0.75 x 10 5 


(4) 




11.7yL/m 


4.40 + 0.44 x 10 5 


(4) 


VLA 4 


6 cm 


0.65 + 0.13 


(1) 




3.6 cm 


0.99 + 0.02 


(2) 


VLA 5 


3.6 cm 


10.9 + 0.22 


(2) 



7 mm, a model image was made of the best-fitting envelope with 
disk mode l fou nd from the SED and image modelling presented 
in Section 15.11 The 7 mm VLA A to D array observations were 
simulated from the model image using the CASA task simdata, 
then combined and imaged. The resultant 7 mm model image did 
not have any significant emission above 3cr = 0.17mJy beam -1 
over the 0.1 arcsec 2 photometry aperture for VLA 3, giving an 
upper limit of 0.57 mJy integrated over the aperture. Hence, we 
estimate that >2.3 mJy of the observed 7 mm emission is from 
ionized gas, and <0.57 mJy is due to dust emi ssion. This is in 
agreement with the results of San na et al. (2012) who found that 
the morphology of the 7 mm emission seen in previous VLA A 
array observations most likely arises from photoionization of the 
outflow cavity walls. 

To measure the spectral index a (where S v oc y a ) of VLA 3, 
which may give further insight into the nature of the emission, 
the 3.6cm and 7 mm data were re-imaged using uv distances 
which were common to both datasets, from 4.5 to 1037 kA, and 
the fluxes remeasured. The integrated fluxes were found to be 
1.43+0.03 mJy for 3.6 cm and 3.3+0. 17 mJy for 7 mm (which 
increased compared to the original 7 mm flux due to a reduc- 
tion in resolution, so that fainter emission was brought above the 
noise level). Therefore the spectral index between 3.6 cm and 
7 mm was found to be 0.5+0.02, where the quoted error in the 
spectral index is due solely to the photometric errors, similar to 
the spectral index expected for an ionized wind, 0.6. However, 
as this value is measured from only two fluxes, it is likely to 
be more uncertain than the photometric errors given. This spec- 




Frequency (GHz) 

Fig. 8. Radio-SEDs of AFGL 2591 -VLA 3 and VLA 1. The SED 
of AFGL 2591 -VLA 3 is show in black error bars, and the re- 
solved fluxes of AFGL 2591 -VLA 3 from lMarengo et al] d2000h 
are shown as red squares. Green squares show the fluxes of 
VLA 1 - errors are not shown as they are smaller than the 
markers. The black and green solid lines show greybody fits 
to the SEDs of AFGL 2591 - VLA 3 between 100 and 85 //in, 
and VLA 1 at mid-IR fluxes, respectively (see Sections [5.2.11 
and 15.2.21) . The green dashed line shows the fit to the long- 
wavelength fluxes for VLA 1 . 



tral index is also an upper limit, as we estimate a fraction of the 
7 mm flux is due to dust emission. By applying the common uv 
distance range above, we derive <0.51 mJy for the contribution 
from dust emission, and therefore >2.8mJy from ionized gas 
emission. Using this lower limit for the ionized gas emission at 
7 mm, we estimate the spectral index to be between 0.4 and 0.5. 

5.2.2. VLA1 

Figure [8] shows the SED of VL A 1 as green squares. Th is source 
was detected in the mid-IR by Mare ngo et al.l (l2000h . as well 
as at several radio frequencies, for which the fluxes are listed 
in Table [7] The spectrum of VLA 1 appears to be flat at radio 
wavelengths, with a fitted spectral index a of 0.0 + 0.03. The fit- 
ted power-law is shown as a dashed line. However, we first need 
to verify whether only a fraction of the flux of VLA 1 is being 
recovered at each wavelength. The flux of VLA 1 measured in 
our 3.6 cm A-D array image is not larger than the fluxes mea- 
sured in the A array only images in previous works (A-D array: 
80 mJy , thi s work; A array only : 82mJy, 94 mJy. pTofani et all 
(1995) and Trinidad et al. (2003) respectively). In fact, there is 
a small decrease in the flux of VLA 1 for the larger uv coverage 
data. Therefore it is likely that most of the flux from VLA 1 has 
been recovered by observations which are sensitive to scales up 
to 7", the largest observable angular scale of VLA A array ob- 
servations at 3.6 cm. This is the case for all of the interferometric 
fluxes shown for VLA 1 , hence the calculated flat spectral index 
of VLA 1 is not likely to be due to instrumental effects. The ob- 
served spectral index is close to that expected from optically thin 
free-free ionized gas emission (a = -0.1). 

The mid-IR emission from VLA 1 is likely to arise from ther- 
mal emission from dust within the HII region. To estimate how 
much VLA 1 contributes to the total flux of the region, a grey- 
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Fig. 9. Continuum emission towards 
AFGL 2591 at 106.7 and 110.5 GHz (-2.8 
and 2.7 mm respectively), shown in both 
contours and greyscale. The rms noise in 
the images is cr =2.2 and 2.1 mJy beam -1 
respectively. The greyscales extend from 
-3 x cr to 65.5 and 54.0 mJy beam -1 . Contours 
are at -3, 3, 4, 5, 7, 10, 15 and 20 x a. The 
crosses show the positions of VLA 1 to 3. The 
synthesised beams are 4.9"x4.1" P.A.= -4.8° 
and 4.3"x 3.5" P.A.= 95° respectively. 



body spectrum with a dust opacity index /? = 2 was fit to the 
fluxes measured by Marengo et al. (2000). The best-fitting grey- 
body has a temperature of 11 1 + 14 K and peaks at 26+3.3 yum 
with a flux of 112Jy, which is uncertain by a factor of two. 
Therefore the thermal dust emission from VLA 1 never con- 
tributes more than -30% to the total flux from the region at 
wavelengths between 8 and lOOyurn and a negligible fraction of 
the observed flux (< 0.5%) at other wavelengths, and therefore 
should not significantl y aff ect the SED and image modelling re- 
sults given in Section 1570 The mid-IR fl uxes for only the re- 
solved central source of AFGL 2591 from Mare ngo et al.l (l2000h 
are plotted on Fig. [8] (red squares) and agree well with the ISO- 
SWS data, which were used to cover the mid-IR to far-IR section 
of the SED. 



5.3. l3 CO, C ls Oand3mm Continuum 

This section details the properties of the dense molecular gas 
surrounding AFGL 2591 inferred from the C 18 molecular line 
and millimetre continuum observations. 

Figure [9] shows the continuum emission from the region at 
106.7 and 110.5 GHz, or 2.81 and 2.7 mm, in both of the wide 
~3 mm continuum bands observed. The morphology of the emis- 
sion is similar to that seen in the mil limetre continuum map s 

JT999), 



at various wavelengths presented by Ivan der Tak et al.. .. 
Ivan der Tak et al.l (120061) . and resolved further by IWang et al.l 
(12012b . We found that the sources in the 3 mm images were 
systematically offset from the 3.6 cm peak positions, which we 
determined to be caused by inaccurate coordinates used for the 
phase calibrator MWC 349 in the mm observations. We therefore 
shifted the continuum and line maps by 0.89" in R.A. and 0.1 8" 
in declinati on to agree with the 7 mm position of MWC 349 A 
reported in lRodriguez et al.l (|2007). In Fig. [9] the peak positions 
of VLA 1 to 3 are shown as crosses. The main two sources that 
are detected are VLA 1 and VLA 3, but there is also extended 
emission towards the position of VLA 2. Although the large un- 
certainties in the two measured 3 mm fluxes do not allow calcula- 
tion of accurate spectral indices for VLA 1 and 2 using these val- 
ues alone, the difference in the spectral index of the two sources 
can be seen from the difference in their relative brightnesses be- 
tween the two images, indicating VLA 1 has a flat spectral index, 
and VLA 3 has a rising spectral index at shorter wavelengths. In 
addition, when we compared both images convolved with the 
beam of the 106.7 GHz image, the morphology of the emission 
was very similar, apart from the clear difference in the flux of 



VLA 3. The sources were fit with 2D Gaussians using the CAS A 
task imfit to find the fluxes, which are given in Table [7] 

Figure [TOl presents C 18 and 13 CO line profiles at the posi- 
tion of AFGL 2591 -VLA 3. The 13 CO observations were miss- 
ing significant flux on extended scales which, in combination 
with self-absorption effects, is likely to be the cause of the dip in 
the 13 CO line profile seen in Fig. [TUJ As there was a significant 
amount of missing flux in the 13 CO line, and the emission was 
distributed incoherently across the maps, it was not possible to 
interpret this emission. Therefore the 13 CO data will not be dis- 
cussed further, however channel maps of the 13 CO emission are 
included in the Online Material as Figs. [TT1 and [T2l 

Figure [13] presents channel maps of the C 18 emission 
towards AFGL 2591 between -9.0 and -2.7 km s" 1 (where - 
5.7 km s" 1 is the rest velocity of the cloud), and Fig. ITU shows 
the C 18 intensity- weighted first moment map. In both Figures 
the positions of VLA 1 through 3 are marked by crosses, how- 
ever in Fig.[l4l the position and size of VLA 1 at 3.6 cm is shown 
instead as a circle. A blue- shifted velocity feature can be seen in 
Fig-Elto the west of VLA 3, and coincident with VLA 2, which 
decreases smoothly in velocity to the west VLA 3, ranging from 
-5.3 to -8.0 km s" 1 . This can be seen in the C 18 channel maps 
(Fig. [13]), which also show that the shape of the emission be- 
comes narrower away from the line centre and towards the west. 
This is most obvious in the -7.7 and -8.0km s" 1 channels. In ad- 
dition, Fig. [14] shows that there is a smaller intensity peak ~6" 
to the south of VLA 3. Here, the velocity instead decreases from 
northwest to southeast from approximately -5 to -6 km s -1 . 

In Fig. Q3] integrated C 18 emission in several velocity 
ranges is compared with the Near-IR Gemini image, and the 
positions of bow shoc ks detected in the K ' band Gemini North 
image (first noted by lPreibisch et al.ll2003l) . The bow shocks are 
shown approximately as white semi-ellipses, as they are too faint 
to be seen directly from Fig. [15] The white dashed line also 
shows the direction of the ionized jet seen in the 3.6 cm images 
(Fig. 0, with RA. 100°. The contours show the high- velocity 
red and blue- shifted gas (red and blue contours respectively). 
Although the "high"-velocity channels referred to in this work 
are high velocity relative to the observed line, it should be noted 
that more extended, higher velocity gas has been detected by 
Hasegawa & Mitchell (1995) using single dish 12 CO observa- 
tions with 14.3" resolution, which extends in velocity from -45 
to 35kms _1 . Therefore, the observations presented here trace 
the higher density but comparatively lower- velocity gas, within 
a few kms" 1 of the line centre. To minimise confusion due to 
overlap in velocity of various components, the contours in Fig. 
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[T5]only show integrated C 18 emission from -4.3 to -3.3 km s 1 
(red) and -8.0 to -7.0km s" 1 (blue). 




-10 -5 

Velocity (km s _1 ) 

Fig. 10. C 18 (black line) and 13 CO (grey line) spectral profiles 
measured at the position of AFGL 2591 - VLA 3 : 20 h 29 m 24 s .88 
+40°11 , 19 ,, .5 (J2000). 



6. Discussion 

6.1. The Jet and Outflow from VLA 3 

In Fig.QJl the extrapolated direction of the ionized jet observed 
at 3.6 cm points directly through the red- shifted C 18 contours 
to the bow- shocks seen in the K' band image. Hence this out- 
lines a coherent picture in which the red- shifted outflow lobe of 
AFGL 2591 consists of a small-scale 4000 AU ionized jet that is 
part of a jet or wind extending out to -0.4 pc where the flow ter- 
minates against the surrounding cloud as bow- shocks. However 
it is interesting to note that the position angle of the jet and 
bow shocks do not exactly align with that of the blue- shifted 
reflection ne bula or the elongation of th e blue- shifted outflow 
observed by lHasegawa & Mitchelll (1 19951) . Therefore other stars 
in the vicinity, such as the powering star(s) of VLA 1 , may be 
causing precession of the jet. If, as well as the emission lying 
along the position angle of the ionized jet, the clump of emission 
at 20 h 29 m 26 s .7 +40° 1 1'23" (J2000) is also included, these three 
clumps may be tracing an arc of emission describing the densest 
parts of the red-shifted outflow lobe, created as the jet precesses. 
If so, the extent of these clumps corresponds to an observed 
opening angle of -40° at their distance from the source, -0.4 pc. 
At the same distance, the observed blue-shifted outflow lobe is 
larger with an opening angle of -60°. The opening angles o f the 
biconical outflow observed by iJimenez-Serra et ail (1201 2|) an d 
the blue shifted outflow cone observed by lSanna et al] ( 20121) . 
on scales <2" or <6700 AU, are -90-110°. The larger observed 
opening angles on smaller scales suggests that the outflow cav- 
ity can not be described by a cone with a constant opening angle, 
but is instead better described by a power-law cavity similar to 
that used in our radiative transfer models. 

Several suggestions regarding the nature of the mm and cm 
continuum emission from VLA 3 have been made by previous 
studies, incl uding a core-halo HII region, an ionized wind with 
a dust disk (Trinidad et al. 2003), or emission from a spheri- 



cal g ravitationally confined HII region (Ivan der Tak & Mentenl 
120051) . At 3.6 cm, the deeper observations presented here clearly 
show that as well as a central compact core, the source also ex- 
hibits a non- spherical jet- like morphology. Therefore, we calcu- 
late several properties of the emission below, assuming it origi- 
nates from a jet. 

Without assuming a specific ionization mechanism , the mass 
loss rate of the jet can be estimated using the model of lReynoldsl 
(1986), which describes the emission from a partially optically 
thick ionized jet (his equation 19): 
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where v is the velocity of the jet, xo is the ionization fraction; 
Li/m p is the mean particle mass per hydrogen atom of the ion- 
ized material, given by 1/(1 + xq); S v is the observed flux at 
the frequency v; a is the spectral index; d is the distance to the 
source; v m is the turn-over frequency below which the emission 
becomes optically thick; 6 is the opening angle of the flow, de- 
fined as the ratio of the projected width to the radius at the base 
of the jet; T is the temperature of the ionized gas; i is the in- 
clination measured from the line-of- sight and F is a function of 
the spectral index an d the dependance of optical depth on radius 
(see lReynoldsl fl986, equation 17). 

Assuming an isothermal, uniformly ionized jet with a den- 
sity gradient of p ~ r~ 2 ; v = 500 km s" 1 for the velocity of the 
j et, found for the b lue- shifted HH objects towards AFGL 2591 
(iPoetzel et al.lll992|); xn = 0.1 (common ly found for low-mass 
sources, e.g. iBacciotti & Eisloffellll999l) : v m = 50 GHz; - 
0.2"/ 1.2", the ratio of the maximum deconvolved width to the 
maximum radius; T = 10 4 K ; a spectral index between 0.4 and 
0.5; and i = 50°, taken from the results of the SEP and near-IR 
image profile modelling presented in Section I5.lt F was found 
to be between 2.3 and 1.8, and the mass loss rate of the jet ob- 
served at 3.6 cm was therefore determined to be in the range 0.77 
- l.OxlO^Moyr" 1 . 

This value is approximately a thousand times larger than the 
jet mass loss rates seen for low-mass protostars, whi ch are com- 
mon ly found to be on the order of 10" 8 M yr _1 (e.g. lPodio et al.l 
2006). The mass loss rate in this jet can also be compared to the 
mass loss rate of the small-scale red- sh ifted molecular flow ob- 
served bv lHasegawa & Mitchelll (1 19951) in 12 CO(J=3-2), which 
corresponds well to the position and direction of the outflow lobe 
suggested by the ionized jet, red-shifted C 18 emission and bow 
shocks discussed in S ection 15.31 For their optically thick case, 
lHasegawa & Mitchelll find a mass loss rate of 6.3 x 10" 5 M©yr _1 
at 1 kpc corresponding to 7.0 x 10" 4 M Q yr _1 at 3.33 kpc. 

Multiplying the ionized jet mass loss rate by the velocity v, 
the momentum rate in the jet is 3.9 - 5.2xl0~ 3 Moyr _1 kms _1 , 
which is v ery similar to the mome ntum rates deter- 
mined by lHasegawa & Mitchelll (1 1995b for both the 
large and small scale red- shifted 12 CO outflows (7.7 and 
7.4xl0- 3 M o yr 1 kms- 1 at 3.33 kpc with i = 50°, scaled from 
7.5 and 7.2xl0 _4 M o yr- 1 kms- 1 at 1 kpc with i = 45°). Thus the 
jet itself would have very close to the required momentum to 
drive the observed larger scale red-shifted outflow. 

The ionization mechanism of the jet can be modelled as a 
plane-parallel shock in a homogeneous neutral 'stellar' wind 
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Fig. 13. C 18 channel map at 0.3 km s" 1 resolution between -9.0 and -2.7 km s" 1 . The rms in the map is <x=0. 1 Jy beam -1 . The map 
peak flux is 1.1 Jy beam -1 . Contours are at -3, 3, 4, 5, 6, 7, 8, 9, 10, 1 1 x <x. The synthesised beam is shown in the bottom left-hand 
corner (4.5"x 3.6", P.A. 93°). A scale size of 50,000 AU is represented by a bar in the bottom right-hand panel. The cross shows the 
position of the central source of AFGL 2591 from the Gemini near-IR J band image. The rest velocity of the cloud is -5.7 km s" 1 . 



(e.g. ICuriel et al.lll989l: lAngladal[T996h . The momentum rate P 
of the jet can be expressed as: 



^M© yr _1 kms 



-l 



3.13 x 10" 



1 



i- 



S v d 2 



\mJy kpc 

j , -0.45 



200 kms -1 
o.i 



0.32 



/ 7 \ ( y \ ( T \ 
\10 4 K/ \5GHz] \l-e- T ) 

(9) 



where rj is the shock efficiency or ionization fraction, found t o 
be ~0.1 for low-mass sources (e.g. Bacciotti & Eisloffel 1999), 
v+ is the initial velocity of the stellar wind or jet, taken to be 
500 kms -1 , and r is the optical depth of the emitting gas. Here, 
the flux S v is measured at 3.6cm. The optical depth can be de- 



termined using equation 6 of lAnglada et al] ([1998): 

In [(1 -e~ Tl )/(l -g- Tl(V1 ^ 2,1 )] 



a = 2 + 



In (vi/v 2 ) 



(10) 



Using the estimated spectral index of 0.4 - 0.5, the optical 
depth of the emission at 3.6 cm was estimated to be between 
2.1 and 2.6, giving a value for the required momentum rate 
in the jet of 2.8 - 3.4xl0 _2 M o yr- 1 kms- 1 . Therefore the re- 
quired momentum rate in the jet is larger than the value of 3.9 
- 5.2xl0 -3 Moyr^kms -1 calculated above by a factor of ~5 to 
9, implying that the 3.6 cm emission cannot be solely caused by 
shocks in the jet. However, if the shock efficiency is increased to 
near unity, the momentum rate of the jet would be sufficient to 
ionize it. As this is unlikely, photoionization by the central star 
therefore must provide a significant fraction of the emission at 
3.6 cm. 
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Fig. 14. Shown in colour-scale: intensity-weighted first moment 
map of the C 18 emission towards AFGL 2591 . The correspond- 
ing velocities are shown on the colour-bar on the right. The black 
open circle represents the position and size of VLA 1 at 3.6cm. 
The crosses indicate the position of VLA 2 and VLA 3. A scale 
size of 20,000 AU is represented by a bar in the upper left-hand 
corner. Shown in contours: the C 18 integrated intensity map 
from -2.3 to -8.7 km s" 1 . The peak flux of the integrated inten- 
sity map is 2.5 Jy beam -1 kms -1 . Contours are at -3, 3, 4, 5, 6, 
7 x <r=0.33 Jy beam -1 kms -1 . The synthesised beam is shown 
in the lower right-hand corner (4.5"x 3.6", P.A. 93°). The rest 
velocity of the cloud is -5.7 km s -1 . 



Turning to the C 18 observations, the blue-shifted emission 
seen in Fig[l4] and the blue contours in Fig. [15] may be trac- 
ing a blue-shifted outflow from VLA 2 or VLA 3. If this is the 
case, the outflow can be seen to be narrower and more colli- 
mated at more negative velocities, seen in Fig.[T3j In addition, a 
"Hubble-like" velocity trend (i.e. an syst ematic increase flow ve- 
locity with distance from the source, e.g. Lada&Fich 1996) can 
be seen across the blue-shifted emission in Fig. [141 Therefore 
there is a high- velocity collimated component at the centre and 
front of this blue- shifted flow, which smoothly transitions into a 
surrounding lower- velocity wide-angle wind. 

Although the integrated C 18 emission peaks on VLA 2, 
thus making it a possibility that VLA 2 is driving the out- 
flow, there are several lines of evidence which suggest this 
is not the case. Firstly, the HII region surrounding VLA 2 
appears to be more evolved than VLA 3, and hence VLA 2 
is less likely to be driving an outflow. Secondly, the emis- 
sion at the position of VLA 3 is at the rest velocity of the 
cloud (5.7 kms -1 ) instead of VLA 2, where it is approximately 
6.3 kms -1 . Thirdly, the blue-shifted emission seen in C 18 
agrees very well with the high- veloci t y larg e-scale emission 
seen in Fig. 3 of Hasegawa& Mitchell (1995). Therefore, if it 
were the case that the blue-shifted C 18 emission was created 
by VLA 2, this source would also be responsible for the large - 
scale blue- shifted outflow seen by Has egawa & Mitchell (1 19951) . 
However, because this source has a much lower total luminosity 
compared to VLA 3, which dominates the emission at all wave- 
lengths, and as this would leave VLA 3 without a corresponding 
blue-shifted outflow to the red-shifted one observed, we consider 
it more likely the blue-shifted emission is the inner part of the 
large-scale blue-shifted outflow of AFGL 2591 -VLA 3. 

Assuming that the observed blue-shifted C 18 emission is 
the inner part of the entrained large-scale blue- shifted outflow of 



Table 8. Derived properties of the inner part of the blue- shifted 
outflow traced by the C 18 O emission 

Parameter Value 
Velocity range (km s _1 ) 
Deprojected outflow length, / (AU) 
Outflow mass, M outflow (M Q ) 
Outflow momentum, M outflow V (M© km s _1 ) 
Outflow energy, 1 /2 M outflow V 2 (M© km 2 s~ 2 ) 
Outflow energy, l/2M outflow V 2 (erg) 
Average velocity relative to v LSR , (V) (km s _1 ) 
Dynamical timescale, t dyn (yr) 



-8 to -5.3 

61,000 

133 

164 

178 

3.54x 10 45 
-1.2 

2.4x 10 5 



AFGL 2591 -VLA 3, the outflow properties were calculated, and 
are given in Table [8] The blue-shifted outflow mass was found 
assuming the gas was optically thin (derived in Appendix A), 



M outflow = 2.75 x 10" 



2 CO 



c 18 o 
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where [^§§] is the abundance ratio of 12 CO to C 18 0, T is the 
temperature of the gas assuming local thermodynamic equilib- 
rium, dkpc is the distance to the source in kpc, and J S v dv is the 
integral of the flux over the velocities of the blue- shifted emis- 
sion. The flux in each channel was found by integrating the emis- 
sion within an irregular aperture around only the blue- shifted 
outflow emission (i.e. not including emission from the resolved 
peak to the south of the central position of AFGL 2591). The 
temperature was assumed to be 25 K fro m the greybody fit to 
the far-IR and sub-mm fluxes in Section 15.2.11 and the abun- 
dance ratio [ 12 co/c 18 o| was calculated using the results given 
in Wilson & Roo 3 {[994), where D GC is the galactocentric dis- 
tance in kpc, 



12, 



CO 



c 18 o 



(58.8 ±11.8)L> + (37.1 ±82.6), 



(12) 



giving a value of = 538 for AFGL 2591 for a galactocen- 

tric distance of 8.5 kpc. 

It should be noted that the shape of the C 18 line profile 
shown in Fig. [TO] may indicate that this line is partially optically 
thick, so that the estimates of the outflow properties depending 
on the derived mass in Table [8] will be lower limits. In addition, 
as the C 18 O emission only traces the densest inner parts of the 
flow, it only represents a small fraction of the mass and momen- 
tum from the entire large-scale outflow. 

To find the momentum M 0Ut fl 0W V and kinetic energy 
\ ^outflow y 2 in the flow, the final term of equation (ITTI) was re- 
placed by J S v (v - v LSR )/ sin i dv or \ J S v ((v - v LSR )/ sin i) 2 dv 
respectively, where v is the velocity of each channel, v LSR =- 
5.7 km s" 1 is the rest velocity of the cloud, and i is the inclination 
of the flow relative to the line of sight, taken to be 50° from the 
results of the SED modelling presented in Section 1570 

The mass weighted average velocity of the flow was cal- 
culated as the momentum over the mass in the flow: (V) = 
outflow V/Af outflow » an d the projected length of the flow was mea- 
sured as the distance from the central source of AFGL 2591 
to the furthest 3 cr contour from the source in the blue- shifted 
emission shown in Fig. Q31 which was found to be ~14.2" or 
47300 AU at 3.33 kpc. Therefore the actual length of the flow 
/ was calculated as / = / pro j / sin i. The dynamical timescale of 
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Fig. 15. Gemini North three-colour JHK' image 
of AFGL 2591 overlaid with contours of C 18 
emission integrated over -4.3 to -3.3 kms -1 
(red) and -8.0 to -7.0 km s" 1 (blue). The three 
arcs in the east of the image show the positions 
of bow shocks visible in the Gemini North im- 
age, which are too faint to be seen directly in 
the Figure. The dashed line shows the direc- 
tion of the ionized jet observed at 3.6 cm (P. A. 
100°), and yellow diamonds show roughly 
the peak positions of th e H 2 knots observed 
in iTamura & Yamashital dl992l) . Red-shifted 
contours are -3,3,4,5,6 x 0.1 Jybeam -1 kms -1 
and blue-shifted contours are -3,3,4,5,6,7 x 
0.1 Jybeam -1 kms -1 . The synthesised beam is 
shown in the lower left-hand corner (4.5" x 
3.6", P.A. 93°). Near-IR image stretch as in 
Fig.0 



the flow was determined using the equation t^ yn = //(V), how- 
ever it should be noted that, as with any outflow, this value may 
not be an exact measure of the time taken to displace the gas, 
as the gas may have been set in motion at some distance from 
the central source. In addition, this assumes a constant outflow 
velocity, but the gas is observed to be accelerating away from 
the source. Nevertheless, this provides a characteristic timescale, 
which can be compared to other studies. Comparing to dynam- 
ical timescales found for large scale flows toward forming mas- 
sive stars (of order 10 4 -10 5 yr for L > 10 4 L^. iBeuther et aD 
2002), the timescale for the inner part of the flow seen in C 18 
is at the upper end of these values. Scaling the blue- shifted dy- 
nam ical timescale of 6.3 x 10 4 y r for the 12 CO outflow observed 
by lHasegawa & Mitchelll (1 19951) to the new distance of 3.33 kpc, 
we obtain a dynamical timescale of 2.1 x 10 5 yr, similar to the 
value we find for the smaller- scale C 18 outflow. 



6.2. Cluster Properties 

The Section below deals with the overall properties of the form- 
ing cluster surrounding AFGL2591-VLA3. The maximum en- 
velope radius for AFGL 2591-VLA 3, found by our modelling in 
Section O is on parsec scales (180,000 AU and 420,000 AU or 
0.87 and 2.0 pc for the envelope with and without disk models 
respectively). This is also corroborated by th e large-scale enve- 
lope observed in the 12 CO and 13 CO maps of Ivan der Wiel et aLl 
(120111) . as well as in the 80" or 270,000 AU diameter d i sk-lik e 
structure seen in the CS(l-O) maps of I Yamashita et al.l (Il987l) . 
This coherent parsec-scale structure, containing at least five 
forming stars (VLA 1 through 5), therefore constitutes a cluster- 
forming clump. 

The spectral types of the stars required to power VLA 1, 2, 
4 and 5 were calculat ed using the 3 . 6 cm fluxes, by the same 
method as outlined in Joh nston et al.l (120091) . assuming they are 
optically thin HII regions with one powering star, giving spectral 
types of B0.5, Bl, B2 and Bl at 3.33 kpc for VLA 1, 2, 4 and 5 
respectively. As VLA 3 has a spectral index of between 0.4 and 
0.5, its emission is not optically thin, and thus such an estimate 



would not be accurate. As the morphology is consistent with a 
compact core plus jet, it may be that some of its emission is not 
from an HC HII region, but is instead created by shock-ionized 
gas. However, assuming photoionization dominates the emission 
from VLA 3, the spectral type of the powering star of VLA 3 was 
estimated to be 08-09 under the assumption of a photoionized 
wind, using the mass loss rate of the jet found in Section 16.11 
0.77 - 1 .Ox 10~ 5 Mpyr" 1 , the assumed wind velocity: 500 kms -1 , 
and Table 4 of iRodriguez & Cantol (Il983l) . For the above spec- 
tral type estimates, we note that the fluxes in Table [6| may only 
represent a fraction of the total flux at these wavelengths, due to 
the finite range in spatial scales which the observations probe, 
and therefore they provide lower limits. 

Following IS anna et~ai1 (120121) . we take the mass of the most 
luminous star VLA 3 to be 38 M Q . Assuming sources VLA 1 
to 5 are cluster members, the lowest luminosity powering star 
(VLA 4 which is B2) corresponds to a ZAMS mass of ~9M , 
found using the ZAMS lumi nosities of iPanagial (1 19731) and the 
theoretical HR diagram of iMaeder & Mevnetl (119891) Thus if 
this cluster has a fully populated IMF, extrapolating a Krupa IMF 
(lKroupal200ll) down to 1 M suggests there should be ~ 100 stars 
in this cluster with mass above 1 M . 

The roughly circular shape of VLA 1 (with a radius of ~2" or 
-6700 AU at 3.33 kpc) and its optically thin spectral index, sup- 
ports the view that it is a UC HII region excited by one or more 
massive stars. However, an exciting star for VLA 1 is not seen in 
either the 2MASS or Gemini North near-IR images. As suming a 
single star, the ZAMS B0. 5 absolute magnitude given in Panagia 
(Il973l) of -2.8, and a distance of 3.33 kpc, the spectral type of 
the calculated exciting star of VLA 1 corresponds to a 2MASS 

band magnitude of 11.1. The detection limit o f the 2MASS 
point source catalog in K s band is 14.3 mag (ISkrutskie et al.l 
120061) . which indicates that there must be more than 3.2 mag of 
extinction at band for this source to be obscured, correspond - 
ing to 26 mag of extinction in the V band (ICardelli et al.ll 19891) . 
These high extinctions, in addition to the indentation that VLA 1 
has made in the ambient cloud traced by near-IR emission, sug- 
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gest that VLA 1 lies within the same cloud core as AFGL 2591- 
VLA3. 

The dent in the cloud at the position of VLA 1 provides sup- 
port to the idea that this UC HII region is at the same distance, 
and is expanding into the material illuminated by the central 
source of AFGL 2591. Therefore it seems likely that the object 
at approximately 20 h 29 m 2 4 s .6 +40° IT 16 " 4 (J2 0QQ), suggested 
to be an edge-on disk by iPreibisch et al.l (120031) . is instead the 
brightened rim of the expanding UC HII region, where it is in- 
teracting with the surrounding cloud. Hence, through its expan- 
sion, VLA 1 may have been able to affect the formation of other 
members of the cluster. 

The secondary peak in Fig. [T4] shows a velocity gradient 
across it in roughly the NW-SE direction. There are no sources 
at this position in any of the images, limiting the possibility 
that this emission is from an outflow of a nearby object. Instead 
the velocity of the emission appears to decrease radially away 
from VLA 1, which is marked as a circle in Fig.[l4j Indeed, Fig. 
|7] indicates that VLA 1 is preferentially expanding towards the 
southeast, into possibly lower density gas. Hence the origin of 
this velocity gradient in Fig. [14] may be due to the expansion of 
the gas around VLA 1 towards the observer, with a projected ve- 
locity of ~1 kms" 1 relative to the ambient cloud. Further, it is 
interesting to note that the lack of C 18 emission at the position 
of VLA 1 again suggests that VLA 1 is at the same distance as 
AFGL 2591-VLA3. The lack of C 18 emission towards VLA 1 
could be due to the fact it has swept-up the molecular cloud ma- 
terial around it, or the molecular gas has been dissociated by the 
high temperatures (~ 10 4 K) in the HII region. 



7. Conclusions 

We have carried out a study of the emission from the luminous 
star-forming region AFGL 2591 from near-IR through to cm 
wavelengths, to determine whether the properties of the source 
can be explained by the same processes thought to govern the 
formation of low-mass stars. 

In Section [U we described a Monte Carlo dust continuum 
radiative transfer dust code to model the continuum absorption, 
emission and scattering through two scaled-up azimuthally sym- 
metric dust geometries, the first consisting of a rotationally flat- 
tened envelope with outflow c aviti es, and the second also in- 
cluding a flared disk. In Section l5Jl our results show that in both 
cases, these models were able to reproduce the main features of 
the SED and 2MASS images of AFGL 2591. 

The observed 3.6 cm images presented in Section [5]2] de- 
tect and more finely resolve the central powering source 
AFGL 2591 -VLA 3, which has a compact core plus collimated 
jet morphology. The east component of the jet extends 4000 AU 
from the central source, with an opening angle of < 10° at this 
radius. This jet is also roughly paralle l to th e outflow seen at 
larger scales by lHasegawa & Mitchell (Il995l) . confirming that 
VLA 3 is the powering source of the outflow from AFGL 2591. 
The multi-configuration 7 mm image detects both VLA 1 and 
VLA 3. At this wavelength, VLA 3 does not show a jet mor- 
phology, but instead presents compact (< 500 AU) emission, 
some of which (<0.57mJy of 2.9 mJy) is estimated to be due 
to dust emission. Comparison with ancillary 0.3-0.4" Gemini 
near-IR images show that VLA 3 is coincident with the illumi- 
nating source of AFGL 2591, at the apex of the observed near-IR 
reflection nebula. Further, VLA 2 is coi ncident with the close bi- 
nary first seen by lPreibisch et al.l (120031) . and VLA 4 and a newly 
detected source VLA 5, are coincident with near-IR sources. 



By fitting greybodies to the SEDs of VLA 1 and VLA 3 at 
wavelengths shorter than 1 mm, the contribution of VLA 1 to 
the total SED from AFGL 2591 was determined to be < 30% 
at wavelengths between 8 and 100/zm, and to be negligible 
at other wavelengths. Therefore contamination by VLA1, or 
fainter sources in the region (e.g. VLA 2), should not signifi- 
cantly change the emergent SED of the central powering object 
which dominates the luminosity from the region, AFGL 2591- 
VLA3. 

The nature of the emission from VLA3 was also investigated. 
The spectral index of VLA 3 was calculated (between 3.6 cm and 
7 mm) to be bet ween 0.4 and 0.5, similar to that of an ionized 
wind. In Section ETTl the 3.6cm emission was modelled as an 
ionized jet, from which the mass loss rate in the jet was deter- 
mined to be in the range 0.77 - l.OxlO -5 M yr _1 . In addition, 
it was found that the jet may have enough momentum to drive 
the larger-scale flow. However, assuming a shock efficiency of 
10%, the momentum rate of the jet was found to not be suffi- 
cient to ionize it via only shocks. Thus a significant portion of 
the 3.6 cm emission is due to photoionization. 

In Section 15.31 the observed C 18 emission towards the 
source was found to be tracing the densest accelerated or en- 
trained material. The blue-shifted emission, to the west of VLA 3 
and coincident with VLA 2, shows an increasing velocity gradi- 
ent to the west, which is more collimated at higher velocities, and 
is coincide nt with the inner part of the blue- shifted outflow ob- 
served by lHasegawa & Mitchelll (IJL995). Thus although we can- 
not rule out that the that this emission is an outflow launched 
by VLA 2, there is substantial evidence that the C 18 emission 
is tracing the densest parts of the entrained gas in the large- 
scale outflow from VLA 3. In addition, the dense red-shifted gas 
traced by C 18 lies along the position angle of the ionized jet 
seen at 3.6 cm, directly between the jet and bow shocks seen in 
the near-IR Gemini image. Therefore it is likely this emission 
also traces a small, dense, segment of the red-shifted flow from 
VLA 3. 

The above results show that many of the properties of 
AFGL 2591 -VLA 3 resemble those of low-mass protostars. For 
example, the geometry of AFGL 2591 determined from mod- 
elling the SED and near-IR images of the source is consistent 
with that of a rotationally flattened envelope with and without 
a flared disk. Furthermore, the outflow of AFGL 2591 -VLA 3 
is comprised of both a collimated ionized jet, and a wide an- 
gle wind, a property s een in the outflo w phenomena of many 
low-mass sources (e.g. lArce et al.ll2QQ7h . However, although the 
emission at 3.6 cm from AFGL 2591 -VLA 3 has the morphology 
of a jet, it is unlikely that it can be explained solely by shocks in 
a neutral wind or outflow. Therefore some part of the compact 
emission from the star may instead be provided by a hypercom- 
pact HII region. Yet if this is the case, its presence has not dis- 
rupted the accretion or outflow processes of AFGL 2591 -VLA 3. 
Thus, in this manner, the formation of the central dominant ob- 
ject VLA 3 of AFGL 2591 does not appear to be significantly 
different to that of low-mass protostars. However, it is important 
to note that this star is not forming in isolation, evidenced by 
the four other cm sources observed within the bounds of the en- 
velope determined for AFGL 2591 -VLA 3, with a radius found 
from our modelling to be on parsec scales. Therefore another 
way to view this picture is that AFGL 2591 -VLA 3 is able to 
source its accreting material from the shared gas reservoir of a 
small cluster while still exhibiting the phenomena expected dur- 
ing the formation of low-mass stars. 
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Appendix A: Mass determination from C 18 
emission 



Starti ng from equation Al from the appendix in IScoville et al.l 
(11986b : 



N total - 



3k 



JiB e W^V)jkT 



T + hB e /3k 



8n 3 B e Li 2 (Ji + 1) (1 - exp (-hv/kT)) 



f 

Jo 



Ty dV 

(A.l) 

where k is the Boltzmann constant, B e is the rotational constant 
for the molecule, \i is the dipole moment of the molecule in e.s.u. 
(1 Debye=10~ 18 e.s.u.), h is the planck constant, // is the rota- 
tional quantum number for the lower energy level, T is the ex- 
citation temperature assumed to characterise the populations in 
all of the energy levels of the gas, v is the frequency of the emis- 
sion and r v is the optical depth over the Doppler line profile as a 
function of velocity v. 

The solution to the equation of radiative transfer (assuming 
background terms are negligible) is 



hv/k 



- Tu ). 



(A.2) 



e hv/kT _ I 

Multiply ing e quation IA.1I by the ratio of the LHS to RHS 
of equation (I A. 21) and simplifying the terms containing e~ hv/kT 
gives 



N total ~ 



N tota i T b (e hvlkT ~ 1) 



hv/k(l 



3k 2 e hBeWi+WT T + hB e /3k 



8n 3 B e Li 2 hv {J i + 1) 



-hv/kT 



r T b Ty 

Jo (1-^ 



- dv. 
(A.3) 



The total mass of gas in the source is given by: 



M, 



H 2 



Ntotal 



H 2 



12 C0 



2 C0 



C 18 



nO 2 



^c m H 2 — d 



(A.4) 



where [it^] is the H 2 to 12 C0 abundance ratio, [^§§] is the 
12 CO to C 18 abundance ratio, ii G is the mean atomic weight 
of the gas, m H is the molecular mass of H 2 , 6 is the angular 
diameter of a uniform disk source, and d is the distance to the 
source. 

Substituting N ma i (equation IA.3b into equation dA.4b we 
find: 



H 2 



3k 2 e hB e Wi+WT T + hB e /3k 



8n 3 B e fi 2 hv (Ji + 1) 

h 2 ir 12 co 



-hv/kT 



Jo 



dv 



12 C0 



c 18 o 



no 2 



(A.5) 



Constants used as well as several values specific to the (1- 
0) transition of C 18 are lis ted below, expressed in c.g.s. units. 
Values for B e are taken from lRosenblum et all ([1958), and values 
for ^ are taken from the JPL line catalog @. 

Ji = 

v= 1.09782182 x 10 11 Hz 

B e = 5.5135449 x 10 10 cycles/s 

ii = 0.11079 Debye = 0.11079 x 10" 18 e.s.u 



H 2 



= 10 4 



12 co 

fi g - 1.36 

m Hl = 3.34524316 x 10" 24 g 



Entering these into equation (IA.5b we obtain: 



M H2 =1.88 x 10" 7 



"12 



CO 



c 18 o 



n2 j2 
u arcsec u kpc 



Jo d-^ 



(T + 0.882) 

e -521/T 

dv M© 



(A.6) 



Similarly, the mass of gas may be obtained from the line flux 
density S v , integrated over a map, 



2k v z 



we find an alternative expression for the mass of the gas: 



(A.7) 



M H2 = 2.75 x 10 -5 



12 
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d 2 
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r SyTy 

Jo (l-^ Ti 



C 18 
) 



(T + 0.882) 

e -521/T 

dv M 



(A.8) 



where S v is measured in Jy. 

If the assumption is made that the emitting gas is optically 
thin (r « 1) over the line, equation (IA.8I) simplifies to 



M H2 = 2.75 x 10" 



12, 



CO 



c 18 o 



(T + 0.882) 2 

p-521/T kpc 



Jo 



S v dv M fl . 



(A.9) 



The abundance ratio [ 12 co/c 18 o l can be calculated using 
the results given in Wilson & Roodl (11994 . where D GC is the 
galactocentric distance in kpc: 



12 



CO 



c 18 o 



(58.8 ± ll.S)D GC + (37.1 ± 82.6). (A.10) 



http://spec.jpl.nasa.gov/ftp/pub/catalog/catdir.html 
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Fig. 11. 13 C0 channel map at 0.3 km s" 1 resolution between -12.0 and -5.7 km s" 1 . The rms in the map is <x=0.1 Jybeam" 1 . The 
map peak flux is 1.3 Jybeam -1 . Contours are at -3, 3, 5, 10, 15 x cr. The synthesised beam is shown in the bottom left-hand corner 
(4.4"x 3.7", P.A. 96°). A scale size of 0.5 pc is represented by a bar in the bottom right-hand panel. The cross shows the position of 
the central source of AFGL 2591 from the Gemini near-IR J band image. The rest velocity of the cloud is -5.7 km s" 1 . 
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Fig. 12. 13 CO channel map at 0.3 kms -1 resolution between -5.7 and 0.6 kms -1 . The rms in the map is <r=0.1 Jy beam -1 . The map 
peak flux is 1.3 Jy beam -1 . Contours are at -3, 3, 5, 10, 15 x cr. The synthesised beam is shown in the bottom left-hand corner (4.4"x 
3.7", P.A. 96°). A scale size of 0.5 pc is represented by a bar in the bottom right-hand panel. The cross shows the position of the 
central source of AFGL 2591 from the Gemini near-IR J band image. The rest velocity of the cloud is -5.7 km s" 1 . 



